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DETAILED ACTION 

1. Claims 1-19 have been examined. 

Drawings 

2. Applicants' formal drawings have been reviewed and approved by the PTO Draftsperson. 

Claim Objections 

3. The following is a partial quotation of 37 CFR § 1.75: 



(i) Where a claim sets forth a plurality of elements or steps, each element or step of the claim should be separated 
by a line indentation. 



3.1 Claims 1-19 are objected to under 37 CFR § 1.75(i) because each element of each 
claim is not separated by a line indentation. 

Claim Rejections - 35 U.S.C. § 112, Second Paragraph 

4. The following is a quotation of the second paragraph of 35 U.S.C. 1 12: 

The specification shall conclude with one or more claims particularly pointing out and distinctly claiming the 
subject matter which the applicant regards as his invention. 

4,1 Claims 12-19 are rejected under 35 U.S.C. 1 12, second paragraph, as being 
indefinite for failing to particularly point out and distinctly claim the subject matter which 
Applicants regard as the invention. 
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4.2 Regarding independent claim 12, the claim preamble is directed to a "method for 
analyzing a design" but the method steps appear directed to mathematical operations for 
approximating a performance surface. Because none of the method steps appear connected to 
analyzing a design, the scope of the claim is unclear. 

4.3 Claims 13-19 are dependent on claim 12 and rejected using the same analysis. 

Claim Rejections - 35 US. C. § 101 

5. The following is a quotation of 35 U.S.C. 101: 

Whoever invents or discovers any new and useful process, machine, manufacture, or composition of matter, or 
any new and useful improvement thereof, may obtain a patent therefor, subject to the conditions and requirements of this 
title. 

5.1 Method claims 7-19 are rejected for reciting a process comprising an abstract 

idea. 

5.2 Regarding independent claim 7, this claim is directed to "a method for analyzing a 
computer generated model," and the steps recited in claim 7 describe the abstract idea of 
performing a mathematical analysis on a variable. Regarding independent claim 12, this claim is 
directed to "a method for analyzing a design," and the steps recited in claim 12 describe the 
abstract idea of performing mathematical operations on a performance surface. 

The steps recited in these claims do not: (1) recite data gathering limitations or post- 
mathematical operations that might independently limit the claims beyond the performance of a 
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mathematical operation; or (2) limit the use of the output to a practical application providing a 
useful, concrete, and tangible result. 

5.3 Dependent claims 8-11 and 13-19 are rejected using the same analysis. 

Claim Rejections - 35 U.S.C §102 
6. The following is a quotation of the appropriate paragraphs of 35 U.S.C. 102 that form the 
basis for the rejections under this section made in this Office action: 

A person shall be entitled to a patent unless - 

(a) the invention was known or used by others in this country/or patented or described in a printed publication in 
this or a foreign country, before the date of invention thereof by the applicant for patent. 

6.1 Claims 1-1 1 are rejected under 35 U.S.C. 102(a) as being anticipated by Ye et al, 
"Algorithmic Construction of Optimal Symmetric Latin Hypercube Design, " Journal of 
Statistical Planning and Inference, Vol. 90 No. 1, pp. 145-159 (July 2000). 

6.2 Regarding claim 7, Ye et al teaches a method for analyzing a computer generated 
model, including: 

receiving the computer generated model [model of diameter of a cyclone, page 22 
equation (3); 

creating at least one variable [seven variables used to predict the diameter of the cyclone; 
see page 22]; and 
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probabilistically analyzing the computer generated model by the use of the at least one 
variable [range of inputs used in modified Latin Hypercube Designs and mean squared error for 
each design, Table 10 and corresponding text at pages 22-25]. 

6.3 Regarding claims 8-9, the method of Ye et al samples a performance surface 
corresponding to 400 randomly selected sites and a calculated mean squared error to analyze the 
design, with the determination of the sampled performance space adequately approximated by 
the lowest mean squared error. See Table 10 and corresponding text. 

6.4 Regarding claims 10-1 1, the method of Ye et al samples a performance surface 
corresponding to mean squared error to analyze the design, with the determination of the 
sampled performance space adequately approximated by the lowest mean squared error. See 
Table 10 and corresponding text. 

Additionally, Ye et al teaches establishing a prediction model for both control and noise 
factors, and then given the distribution of noise variables, estimating the variation of the response 
for each combination of control variables using the model. See page 25 paragraph 1. 

6.5 Regarding system claims 1-4, these claims correspond to method claims 7-9 and 
are anticipated using the same analysis. 

6.6 Regarding system claims 5-6, these claims correspond to method claims 10-11 
and are anticipated using the same analysis. 
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Conclusion 

7. The prior art made of record and not relied upon is considered pertinent to Applicants' 
disclosure. Reference to Ali et al, U. S. Patent 6,405,344 issued 1 1 June 2002 and filed 14 May 
1999, is cited as teaching a method for performing design trade-off based on a total score. 

Reference to Ostrowski et al, U. S. Patent 6,377,908 issued 23 April 2002 and filed 14 
May 1999, is cited as teaching a method for optimizing transfer function outputs in which inputs 
are perturbed. 

Reference to Bailey et al, "Using Response Surfaces to Improve the Search for 
Satisfactory Behavior in System Dynamics Models," System Dynamics Review, Vol. 16 No. 2, 
pp. 75-90 (2000), is cited as teaching an augmented robust concept exploration method. 

Reference to Hamada, "Using Statistically Designed Experiments to Improve Reliability 
and to Achieve Robust Reliability," IEEE Transactions on Reliability, Vol. 44 No. 2, pp. 206- 
215 (June 1995), is cited as teaching Taguchi experimental designs. 

Reference to Boning et al, "DOE/Opt: A System for Design of Experiments, Response 
Surface Modeling, and Optimization Using Process and Device Simulation," IEEE Transactions 
on Semiconductor Manufacturing, Vol. 7 No. 2, pp. 233-244 (May 1994), is cited as teaching the 
DOE/Opt system that integrates design of experiments, response surface model generation, and 
nonlinear constrained optimization. 

8. This Office action has an attached requirement for information under 37 CFR 1.105. A 
complete reply to this Office action must include a complete reply to the attached requirement 
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for information. The time period for reply to the attached requirement coincides with the time 
period for reply to this Office action. 

9. Any inquiry concerning this communication or earlier communications from the 
Examiner should be directed to Samuel Broda, whose telephone number is (703) 305-1026. The 
Examiner can normally be reached on Mondays through Fridays from 8:00 AM - 4:30 PM. 

If attempts to reach the Examiner by telephone are unsuccessful, the Examiner's 
supervisor, Kevin Teska, can be reached at (703) 305-9704. The fax phone number for the 
organization where this application or proceeding is assigned is (703) 872-9306. 

Any inquiry of a general nature or relating to the status of this application or proceeding 
should be directed to the group receptionist, whose telephone number is (703) 305-3900. 

SAMUEL BRODA, ESQ. 
PRIMARY EXAMINER 
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REQUIREMENT FOR INFORMATION UNDER 37 C.F.R. § 1.105 
ATTACHMENT TO PAPER NO. 200408 

1. Applicants and the assignee of this application are required under 37 CFR 1.105 to 
provide the following information that the examiner has determined is reasonably necessary to 
the examination of this application. 

2. In response to this requirement, please provide copies of each publication (published 
prior to the filing date) which any of the Applicants authored or co-authored and which describe 
the disclosed subject matter of: 

(1) Latin hypercube designs; or 

(2) entropy analysis. 

3. In responding to those requirements that require copies of documents, where the 
document is a bound text or a single article over 50 pages, the requirement may be met by 
providing copies of those pages that provide the particular subject matter indicated in the 
requirement, or where such subject matter is not indicated, the subject matter found in 
Applicants' disclosure. 

4. The fee and certification requirements of 37 C.F.R. § 1 .97 are waived for those 
documents submitted in reply to this requirement. This waiver extends only to those documents 
within the scope of this requirement under 37 C.F.R. § 1.105 that are included in the Applicants' 
first complete communication responding to this requirement. Any supplemental replies 
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subsequent to the first communication responding to this requirement and any information 
disclosures beyond the scope of this requirement under 37 C.F.R. § 1.105 are subject to the fee 
and certification requirements of 37 C.F.R. § 1.97. 

5, The applicant is reminded that the reply to this requirement must be made with candor 
and good faith under 37 CFR 1.56. Where the applicant does not have or cannot readily obtain 
an item of required information, a statement that the item is unknown or cannot be readily 
obtained will be accepted as a complete response to the requirement for that item. 

6. This requirement is an attachment of the enclosed Office action. A complete response to 
the enclosed Office action must include a complete response to this requirement. The time 
period for reply to this requirement coincides with the time period for reply to the enclosed 
Office action, which is THREE months. 
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Algorithmic construction of optimal symmetric 
Latin hypercube designs 



Kenny Q. Ye William Li Agus Sudjianto 

Department of Applied Department of Operations Reliability Methods Department 
Mathematics and Statistics and Management Science 

SUNY at Stony Brook University of Minnesota Ford Motor Company 

Stony Brook, NY 10021 Minneapolis, MN 55455 Dearborn, MI 48121 

Abstract 

We propose symmetric Latin hypercubes for designs of computer exepriment. 
The goal is to offer a compromise bet ween computing effort and design opti- 
mally. The proposed class of designs has some advantages over the regular 
Latin hypercube design with respect to criteria such as entropy and the min- 
imum intersite distance. An exchange algorithm is proposed for constructing 
optimal symmetric Latin hypercube designs. This algorithm is compared 

1 



with two existing algorithms by Park (1994) and Morris and Mitchell (1995). 
Some examples, including a real case study in the automotive industry, are 
used to illustrate the performance of the new designs and the algorithms. 
Key Words: Computer experiment; Maximum Entropy Design; Maximin 
design 

1 Introduction 

One of our recent projects is concerned with the thermal analysis of multi- 
layer electrical traces at a major automotive company. As more and more 
electronic devices are installed in vehicles, the peak temperature of electrical 
traces becomes a major concern in designing the instrument panels. The 
temperature of an electrical trace is largely determined by its width, its 
passing current, strength, and its position in a stack of traces. The goal of 
this project is to provide guidelines for design engineers for width and passing 
current strength of multi-layer electrical traces. Physical experiments are 
inevitably very expensive and time consuming since a set of electrical traces 
has to be assembled in certain configurations for each test and measuring 
the temperature of each trace is difficult. Therefore, finite element analysis 



(FEA) models have been developed to simulate the thermal dynamics of 
electrical traces. 

Using the computer model, the study starts from a simple case, in which 
there are two layers with three traces on each layer. One primal}' interest 
is the interaction between a center trace and an edge trace on two different 
layers since the heat coming off the center trace spreads out and affects the 
temperature at the edge. A center trace on layer* 1 and an edge trace on layer 
2 are selected in the study. The goal is to predict their peak temperatures 
(ill and y 2 ) based on four predictors: the width of the center trace (x { ).. the 
applied current of the center trace (#2), the width of the edge trace (£3), 
and the applied current of the edge trace (j.4). Given a set of ;r,-values. the 
computer model generates a deterministic peak temperature for each trace. 

Though computer experiments are much cheaper and faster than physical 
experiments, each run is still time consuming and expensive. Thus, only a 
small number of combinations of the can be tested. In this case, the exper- 
iment is to be conducted by another company which specializes in thermal 
dynamic computer models, and the budget only allows for 25 runs. A feasible 
approach is to establish a statistical model from the results of the 25 runs 
and then use it to predict peak temperatures for any given combinations of 



the Xi. An optimal Latin hypercube design was chosen for this experiment. 

A Latin hypercube design (LHD) is an n x / matrix in which each col- 
umn is a random permutation of {1....,??} which can be mapped onto the 
actual range of the variables. It has good projection properties on any single 
dimension. Latin hypercube designs have been applied in many computer 
experiments since they were proposed by McKay et ai (1979). In practice, 
a LHD can be randomly generated, but a randomly selected LHD may have 
bad properties and act poorly in estimation and prediction. Another ap- 
proach is to use optimal designs according to some criteria such as entropy 
(Shewn- and Wynn 1987). Integrated Mean Squared Error (IMSE) (Sacks 
et al. 1989), and minimum intersite distance (Johnson et ai 1990). These 
designs have been shown to be efficient for certain models. However, the 
computational cost of obtaining these designs is high. In an attempt to offer 
a compromise between good projective properties of LHDs and a criterion, 
Park (1994) and Morris and Mitchell (1995) proposed optimal Latin hyper- 
cube designs. For an excellent review of design and analysis of computer 
experiments, see Koehler and Owen ( 1996). 

One of the criteria considered in this paper is the entropy criterion, 
first proposed by Shewn' and Wynn (1987) and then adopted by Currin 



et al (1991). The response of a computer model is modeled by F(x) = 
+ ^( x )* ^ which Z(x) is a Gaussian process with mean zero 
and covariance 

/?(s. t) = a 2 exp |-0 £j | Sj - t/| ? 0 < q < 2 (1) 

between two /-dimensional inputs s and t. The entropy criterion is equiva- 
lent to the minimization of -log|/?|, where R is the covariance matrix of the 
design. The parameters 9 and q determine the properties of Z(x). Through- 
out this paper, we set q = 2 so that the correlation between two sites is a 
function of their L 2 distance. 

The construction of an optimal LHD can still be time consuming. For 
example, to generate a maximum entropy 25 x 4 LHD using a columnwise- 
pairwise (CP) algorithm (discussed in Section 3). the whole procedure takes 
3.3 hours on a Sun SPARC' 20 workstation, which appears to be quite long- 
as the size of the design is moderate. The search for a larger design would 
take even longer, and may be computationally prohibitive. This situation 
motivated us to look for alternatives that require less computing time. The 
easiest method is to generate a large number of random LHDs and then 
choose the best, one according to the criterion. For example, the generation 
of 1000 random LHDs takes only 14.7 seconds on the same machine. However, 



the best design obtained from these random designs is usually significantly 
inferior to that produced by the algorithmic search. In our example, the 
entropy value at 0 = 2 of the former is 25.26, compared with the latters 
20.48. To reduce the searching time and still generate competitive designs, 
our approach is to restrict the search within a subset of the general LHD. If 
this subset of designs has some desirable properties with respect to a criterion, 
then selecting a design from this group of designs may be more efficient. 

In Section 2 7 we introduce a new class of LHD, the symmetric Latin hyper- 
cube design, whose geometric property enables us to find optimal LHDs more 
efficiently. Section 3 considers a simple exchange algorithm for constructing 
optimal symmetric LHDs. Its performance is compared with the existing 
algorithmic approaches of Park (199-1). and Morris and Mitchell (1995). Sec- 
tion 4 demonstrates the performance of the new design with an example. A 
summary is given in Section 5. 

2 Symmetric Latin hypercubes 

Our goal is to find a special type of LHD that has some good "built-in* 1 
properties with respect to the optimality of a design. In our definition, a 



C 



LHD is called a Symmetric Latin hypercube design (SLHD) if it has the 
following property: in an n x / LHD with levels from 1 to r?. if (r/i. a. 2 , • — , (ii) 
is one of the rows, then the vector (n + 1 — n + 1 — r/ 2 , *•-,-». + 1 — (h) 
must be another row in the design matrix. In other words, if t,; is a design 
point in a SLHD, then there exists another point tj in the design that is the 
reflection of t 7 through the center. An example of a 10 x 5 SLHD is given in 
Table 1, in which the i th row is the symmetric point of the (n + 1 — i) fh row. 



1 2 3 4 5 

1 6 C 5 9 

2 2 3 2 4 

3 19 7 5 

4 3 4 10 3 

5 7 18 10 

6 4 10 3 1 

7 8 7 1 8 

8 10 2 4 C 

9 9 8 9 7 

10 5 5 6 2 



Table 1: A 10 x 5 symmetric Latin hypercube design 



The symmetry of a SLHD provides some orthogonal properties. That 
is, the estimation of the linear effect of each variable in a SLHD is uncorre- 
lated with all quadratic effects and bidinear interactions. The proof of such 
properties can be found in Ye (1998), in which Orthogonal Latin Hypercube 



Designs (OLHD) are constructed and proposed. OLHDs have the same sym- 
metric properties but also process additional orthogonality which insures the 
zero correlation among estimation of linear effects. Therefore, one can view 
the SLHD as a generalization of the OLHD. However, the number of runs 
in an OLHD has to be a power of two, which increases dramatically as the 
number of factors increases. In the cases that an appropriate OLHD can 
not be found under the constraint of run size, one can consider using SL-HDs 
which have the flexibility of the run size, yet retain some of the orthogonality 
of an OLHD. 

Are SLHDs better than regular LHDs with respect to design criteria? 
Many optimal LHDs reported by Park (1994) and Morris and Mitchell (1995) 
have some symmetric properties. In particular, Morris and Mitchell (1995) 
noticed that a large number of the optimal LHDs they obtained are SLHDs 
and referred them as "fold over designs". This is the first time the SLHD 
is mentioned in the literature. They called for a thorough investigation on 
this phenomenon. Intuitively, the optimal designs are considered to have- 
good space filling properties, and a good space filling design probably lias 
some degree of symmetry. To verify this, we undertook a simulation study to 
compare random SLHDs to random LHDs with respect to both entropy and 
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minimum intersite distance criteria. Table 2 compares the best design among 
the 1000 random SLHDs of size 25 with that of the 1000 random regular 
LHDs. The former has a smaller entropy criterion value of 23.60 at 0 — 2, 
compared with the latter's 25.26. In the table, we also list the minimum 
distances of three designs, a criterion first proposed by Johnson et al (1990) 
and then used by Morris and Mitchell (1995) in constructing optimal LHDs. 
For a design S. the minimum distance d*(S) = min s ^ 6(S - d(s. t), where s 
and t are two design points [i.e., two rows in the design). Both the L[ 
(rectangular) distance L^s^t) — Y?j=i \ s j — t/l an< l ^2 (Euclidean) distance 
L* 2 (s,t) — [Ylj=i(Sj — tj) 2 ]5 of three designs are listed in Table 2. A design 
Si is said to be better than design S 2 if d*(S\) > d*(S 2 ). The number 
of pairs separated by this distance, denoted J, is shown in parentheses in 
the table. If two designs have the same </* values, then the design with 
smaller J value is better. Throughout this paper, entropy and distances are 
computed after LHDs are scaled into [0, 1]'. The levels of scaled LHDs are 
{0^ 7~t- ■ • ••: 1}? which are the same used by Morris and Mitchell (1995). 
Note that Park (1994) used a different scale, {£, f , • • ■ : . 1}. 

Tables 3 and 4 provide more comparisons between regular LHDs and 
SLHDs. We study Latin hypercubes with six different sizes, 25 x 4, 20 x 6. 
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design 


entropy 


min. Ly distance 


min. L-2 distance 


best, random LHD 


25.26 


0.50 (T) 


0.29 (1) 


best, random SLHD 


23.60 


0.63 (2) 


0.35 (2) 


optimal LHD 


20.48 


0.75(1) 


0.50(1) 



Table 2: Comparison of three 25x4 LHDs by (i) entropy criterion (0 = 2); fii) 
smallest L { distance and the number of pairs separated by that distance (in 
parentheses); and (iii) smallest L 2 distance and the number of pairs separated 
by that distance. The "best" is in terms of the entropy criterion. 



LHD SLHD 



Size 


Mean 


Max 


Mean 


Max 


25 x 4 


0.3478 


0.5417 


0.3944 


0.625 


20 x 6 


0.8205 


1.211 


0.8968 


1.316 


50 x 10 


1.316 


1.735 


1.407 


1.816 


200 x 20 


2.9457 


3.4925 


3.0884 


3.5678 


500 x 30 


4.8737 


5.3567 


5.0240 


5.4188 


1000 x 50 


9.4014 


9.9339 


9.6016 


10.2002 



Table 3: Minimum L Y distances of random SLHDs and random LHDs 

50 x 10 , 200 x 20, 500 x 30, 1000 x 50. The sample sizes are 1000 for the 
first three designs and 100 for the last, two designs. SLHDs are consistently 
superior to the corresponding LHDs with respect to both L Y and L 2 distances. 
We also compare the means of minimum distance of LHDs and SLHDs using 
/-tests. The p-values are all smaller than 0.000.1. Therefore. SLHDs are 
statistically significantly better than LHDs with respect to the minimum 
distance criteria. 
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LHP SLHD 

Size Mean Max Mean Max 

25 x 4 0.1943 0.3200 0.2230 0.3727 

20 x 6 0.3875 0.5642 0.4270 0.G31C 

50 x 10 0.51G3 0.7061 0.5492 0.7210 

200 x 20 0.8492 1.0568 0.8879 1.0281 

500 x 30 1.1383 1.2812 1.1812 1.2960 

1000 x 50 1.7327 1.8328 1.7658 1.8689 

Table 4: Minimum L 2 distances of random SLHDs and random LHDs 

These simulation studies have shown the advantages of "picking the win- 
ner" from SLHDs instead of regular LHDs. However, the best. SLHD obtained 
by the "picking the winner' approach is usually inferior to the corresponding 
optimal design obtained by a searching algorithm, as shown in Table 2. In 
the next section, a simple exchange algorithm is presented to search optimal 
SLHDs. 

3 An algorithm and examples 

Tn this section, we review the two existing algorithms proposed by Park 
(1994) and Morris and Mitchell (1995), respectively for constructing opti- 
mal LHDs. Then a eolumnwise-pairwise exchange algorithm (CP) is intro- 
duced in the context, of the construction of optimal SLHDs. The similarities 
and differences between the CP and the other two algorithms are discussed. 
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Through examples, we compare 

(1) the performance of the CP and other two algorithms; 

(2) the optimal SLHDs and the optimal regular LHDs with respect to the 
design criteria and the searching time. 

3.1 Existing algorithms 

To construct optimal LHDs : Park (1994) presented an approach based on 
the the exchanges of several pairs of the elements in two rows. His algorithm 
first selects some active pairs which minimize the objective criterion value 
by excluding that pair from the design. Then for each chosen pair of two 
rows and l 2 , the algorithm considers all the possible exchanges x iin <-> 
z-hji - • ■ ■ • J'nn ^ *i,> Jk . for A' < I and finds the best, exchange among them. 

Morris and Mitchell(1995) adopted a simulated annealing algorithm to 
search for optimal LHDs. They also defined a maximin distance criterion. 
For a given design, define a distance list {d u rf 2 , . . . , //,„.).</•. < rf 2 < . . . < d m . 
in which the are the distinct, values of intersite distances. Let J t be the 
number of pairs of sites in the design separated by d t . Then a design is a 
maxirnin distance design if and only if 

(la) d[ is maximized, and among the designs for which this is true, 
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(lb) J[ is minimized, and among the designs for which this is true, 
(2a) d<2 is maximized, and among the designs for winch this is true, 
(2b) J-2 is minimized, 

and so forth. Morris and Mitchell (1995) also pointed out that, although this 
extended definition is intuitively appealing, it would be better to use a scalar- 
valued criterion as the driving criterion. For this purpose, they proposed a 
familv of functions 



4> P 



Up 



(2) 



where p is a positive integer. Normal]}-, different p values are tried to obtain 
a maximin distance LHD. 

In Morris and Mitchell's algorithm, a search begins with a randomly 
chosen LHD, and proceeds through examination of a sequence of designs, 
each generated as a perturbation of the preceding one. A perturbation D try 
of a design D is generated by interchanging two randomly chosen elements 
within a randomly chosen column in D. The perturbation D trv replaces D 
if it leads to an improvement. Otherwise, it will replace D with probability 
- = exp{ — [o(D try ) — &(D)]/t}, where t is a preset parameter known as the 
"temperature". 
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3.2 Our algorithm 



Li and Wu (1997) considered a class of columnwise-pairwise algorithms in the 
context of the construction of optimal supersaturated designs. A columnwise 
algorithm makes exchanges on the columns in a design and can be partic- 
ularly useful for designs that have structure requirements on the columns. 
Note that each column in an n-run LHD is a permutation of {1, ... ,77-}. At 
each step, another permutation of (1, . . . . n} is chosen to replace a column 
so that the Latin hypercube structure is retained. Therefore, we adopt the 
colummvise-painvise idea in searching for optimal LHDs. However, one im- 
portant change has to be made to accommodate the special structures of 
the SLHD. For a SLHD two simultaneous pair exchanges are made in each 
column to retain the symmetry. For example, suppose a column in a G-row 
SLHD is (L 2, 3, 4, 5, 6)'. If element 1 is exchanged with /. element 6 must be 
exchanged with n + 1 - / (i.e. 7 - /) to keep the design symmetric. The only 
exception is when element. / is exchanged with element /? + 1 — /, which does 
not require a second exchange. The exchange procedure for a SLHD with an 
odd number of rows is slightly different. The center point of the design does 
not participate in the exchange. For example, if a column in a 7- row SLHD 
is ( 1, 2, 3, 4. 5. G, 7)', then element 4 may not to be exchanged with any other 
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element. 

The algorithm for searching optimal SLHD is summarized as follows: 

1. Start with a random SLHD. 

2. Each iteration has / steps. At the /tli step, the best two simultaneous 
exchanges within column / are found. The design matrix is updated 
accordingly. 

3. If the resulting design is better with respect to the criterion, repeat 
Step 2. Otherwise, it is considered to be an "optimal design' 1 , and the 
search is terminated. 

The resulting optimal designs depend largely on the starting designs used in 
the algorithm. Hence, one should repeat the algorithm with several differ- 
ent random starting designs. The best, design among the generated optimal 
designs is chosen to be the final design. 

3.3 Examples 

Example 1 CP vs. Simulated Annealing 

The simulated annealing algorithm proposed by Morris and Mitchell (1995) 
aims at constructing optimal regular LHDs. We modify their algorithm to 



search for SLHDs. Similarly, the CP algorithm discussed in the previous 
section is modified to construct optimal regular LHDs. Both algorithms are 
colummvise-painvise procedures. The simulated annealing algorithm oper- 
ates on a (randomly chosen) column and then considers a (randomly chosen) 
pair in each column. Our proposed CP algorithm resembles the former with 
a very low starting temperature (so that, switches to inferior designs are never 
made). An important difference is that the simulated annealing algorithm 
perturbs the design in a random manner, and our CP algorithm perturbs the 
design in a deterministic manner. 

To compare their performances, we use both algorithms to construct opti- 
mal regular LHDs and optimal SLHDs. Two examples are considered. Table 
5 lists the 12 x 2 maximin distance LHDs and SLHDs generated by both 
algorithms- The simulated annealing algorithm uses 10 starting designs and 
the CP uses 100 starting designs. The driving criterion is 0 p with p = 50. 
To compare the efficiency of the algorithms, the number of searched LHSs is 
also recorded. Bo tit algorithms obtain equally good optimal designs. But the 
CP algorithm searches far fewer LHDs than the simulated annealing algo- 
rithm does. When both algorithms are used to construct 25 x 4 designs, the 
simulated annealing algorithm produces bettor optimal designs than the CP. 
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design 


algorithm 


min. 


dist. 


# of searched LHDs 


LHD 


Simulated Annealing 


.4545 


(16) 


269520 


LHD 


CP 


.4545 


(16) 


44220 


SLHD 


Simulated Annealing 


.4545 


(16) 


2404 1C 


SLHD 


CP 


.4545 


(16) 


14652 



Table 5: Comparison of optimal 12x2 LHDs and SLHDs using two al- 
gorithms, CP (100 starting designs) and simulated annealing (10 starting 
designs). The search criterion is <j> p with p = 50 and L 2 distance. Note: 
using the simulated annealing search in LHD, only two out of 10 starting 
designs result in 0.4545 (1C), compared to all 10 in the SLHD case. 



design 


algorithm 


min. Li dist. 


# of searched LHDs 


LHD 


Simulated Annealing 


.9177 (19) 


1537GC3 


LHD 


CP 


.8750 (6) 


2241900 


SLHD 


Simulated Annealing 


.9583 (36) 


1426985 


SLHD 


CP 


.9177 (6) 


546480 



Table G: Comparison of optimal 25 x 4 LHDs and SLHDs using two al- 
gorithms, CP (100 starting designs) and simulated annealing (10 starting 
designs). The search criterion is <b p with p = 50 and L y distance. 

Therefore, we ma^- conclude that, the systematic search algorithm is better 
for small designs and the simulated annealing algorithm is bet ter for larger 
designs. 

Example 2 CP vs. Park 

Park's algorithm (1994) cannot be easily modified to accommodate the 



property of symmetry. Thus, its comparison with the CP is done through 



construction of the optimal 9x2 regular LHDs, which is discussed in detail by 
Park(1994) to illustrate his exchange algorithms. The CP algorithm is also 
used to construct 9x2 SLHDs. Table 7 compares the optimal designs gen- 
erated by two algorithms with respect to the entropy criterion (0 = L 5, 25), 
along with the total number of searched LHDs for each algorithm. Three 
interesting observations are apparent in this example: 

1. The CP seems to consistently produce better LHDs than Park's al- 
gorithm with respect to entropy. The former also reaches the final design 
slightly earlier since it searches fewer LHDs. In fact, exhaustive searches re- 
veal that the CP produces the global optimum for each value of 9 = 1, 5, 25. 
Our study of constructing LHDs of different sizes shows the same patterns. 

2. Comparisons between optimal LHDs and the corresponding SLHDs 
show that the former are slightly better than the latter but take approxi- 
mately 4 times as much time to search. However, for such a small design, 
it takes so little time (G to 7 seconds on a Sun Sparc 20 workstation) for an 
exhaustive search in the whole LHD class. Therefore, there is no need to 
restrict the search within SLHD class. 

3. At 0 = 25. the global optimal LHD is symmetric. Moreover, it is also 
an orthogonal Latin hypercube as constructed algebraically by Ye (1998). 
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9 


design 


algorithm 


optimal 


average 


# of searched LHDs 


1 


SLHD 


CP 


20.38 


22.16 


5488 




LHD 


CP 


19.1C 


19.29 


20592 




LHD 


Park 


20.01 


20.79 


24132 


5 


SLHD 


CP 


3.09 


4.35 


4768 




LHD 


CP 


2.95 


3.0G 


20052 




LHD 


Park 


3.42 


3.79 


23970 


25 


SLHD 


CP 


0.49 xl0~- 


0.83 xlO -2 


4784 




LHD 


CP 


0.49 x 10 -2 


1.11 xlO" 2 


19080 




LHD 


Park 


1.03 xlO -2 


4.67 xlO- 2 


23580 



Table 7: Comparison of three algorithms for generating 9x2 optimal LHDs 
with 100 random starting designs 

Example 3 LHD vs. SLHD 

We now revisit the case study at the beginning of this paper: the con- 
struction of a 25 x 4 LHD for the thermal analysis of electrical traces. A 
primary motivation of using the SLHD is to reduce the searching time. Since 
the number of possible exchanges of each column in a SLHD is much less 
than that, for a regular LHD, it is expected that the exchange algorithm for 
the SLHD will use much less CPU time. This is confirmed when the algo- 
rithm is applied to the construction of the optimal 25 x 4 SLHD. Without 
the restriction to SLHD, it takes 13.3 hours and 10. G hours for the algorithm 
to terminate using the entropy and o p criteria, respectively. Using the same 
number of starting designs (100), the optimal SLHDs are found only after LG 
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hours and 1.3 hours. The results are summarized in Table 8. Theoretically, 
the global optimal SLHD cannot be better than the global optimal LHD 
since the SLHD is a subset of the LHD. It is seen in our Example 2 that the 
obtained optimal LHDs are globally optimal verified by exhaustive search, 
but they do not always have the symmetrical structure. Morris and Mitchell 
(1995) use an exhaustive search to find maximin distance LHDs for many 
small designs. Not all those global optimal designs are symmetric. In prac- 
tice, a globally optimal design is rarely obtained when the exhaustive search is 
not feasible. In our case, with much less searching time, the optimal SLHDs 
found are actually better than the two optimal LHDs obtained previously 
with respect to both entropy and minimum distance criteria. The maximum 
entropy SLHD has the criterion value of 18.53 compared with 20.48 for the 
previously obtained optimal LHD. The former also has the better (i.e. larger) 
minimum distance of 0.83, with six pairs separated by this distance. Using 
o p as the driving criterion, a SLHD was obtained with four pairs separated 
by the minimum intersite distance of 0.92. which is considerably bettor than 
the optimal maximin distance LHD previously found (six pairs separated by 
0.83). 

Now we revisit Table G and focus on the difference between SLHD and 
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design 


criterion 


entropy 


min. Li dist. 


CPU time(hrs) 


LHD 


entropy 


20.48 


.75 (1) 


13.34 


LHD 


<h 


23.52 


.83 (6) 


10.C3 


SLHD 


entropy 


18.53 


.83 (6) 


1.6 


SLHD 


<h 


19.48 


.92 (4) 


1.28 



Table 8: Comparison of optimal 25 x 4 LHDs vs. SLHDs with respect to 
entropy and cb p . The entropy is calculated for (9 — 2. 

regular LHD. For the simulated annealing algorithm, restricting the search 
within the SLHD did not save much searching time, but the obtained design 
is significantly better with respect to the minimum distance criterion. For the 
CP algorithm, there is a dramatic reduction in searching time after we restrict 
the search within the SLHDs, yet the obtained design is much better. It is also 
interesting to observe that in far less time the CP found a better design within 
the SLHDs than the simulated annealing algorithm found within general 
LHDs. 

4 A robust design simulation example 

One of the goals in computer experiments is prediction. The Kriging method 
was developed in geostatistics and brought into computer experiment by 
Sacks el al (1.989) and Currin et ai (1991) to predict untested sites in the 
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experimental regions. It models the response as a Gaussian process. Given 
a correlation function of the process, the best linear unbiased predictor of y 
at site x is given by 

y(x) =/i + r T fl- 1 (y-/il) ; 

where r is the vector of correlations between x and the design sites x,;. R is 
the correlation matrix among design sites, and y is the vector of the observed 
responses. In this section, an example is used to illustrate the advantages of 
using optimal SLHDs for Kriging methods. 

The example presented here is taken from Mori (1985), which was orig- 
inally presented as a robust design case study. It was later used by Li and 
Wu (1999) to illustrate an integrated approach to parameter design and 
tolerance design. .The original study is concerned with the design of cy- 
clones, which are used to separate solid mass and gaseous mass in chemical 
engineering. There are seven variables whose original values are given by 
(Xi.x^x^x^x^xc^xr) = (0.1,0.3,0.1,0.11.5,10.0,0.75). The relation be- 
tween the response, the diameter of a cyclone (y), and these seven variables 



is 



, /=l74 .42(%— E3_)*«x 



1 - 2.62(1 - 0.36(^)-°-^yV->(2±)i-tc 



xaxr 

The range of the input variables in this example is taken to be 1 ± 10% of 



the original value (see Table 9). 



Unsealed Input Lower limit Upper limit 



x { 0.09 0.11 

x-2 0.27 0.33 

j: 3 0.09 0.11 

x 4 0.09 0.11 

x 5 1-35 1.65 

*c 14.4 17.6 

x 7 0.675 0.825 



Table 9: Range of the inputs in equation (3) 

Experiments with 16 runs are performed using (1) 100 random LHDs; (2) 
100 random SLHDs; (3) maximum entropy SLHDs generated by the CP with 
0 = 0.05,0.5, 1; and (4) a maximin I 2 distance SLHD generated by the CP 
with In each experiment, using the Kriging method with the correlation 
function given in equation (1) with 9 = 0.05,0.1.0.5. we predict Y(x) at 
the same 400 randomly selected sites. The mean squared error (MSE) of 
predictions at these 400 sites was calculated for each experiment. The results 
are summarized in Table 10. First, it can be seen that. r,he MSE is sensitive 
to the 9 used in the Kriging model but is insensitive to the optimal design 
criterion. Second, all of the optimal designs are better than the random 
designs. Third, in this case, SLHDs do not always outperform LHDs. 

In practice, the choice of correlation function in the Kriging model is 



Correlation parameter in kriging model 



Design 


0 = 0.05 


0 = 0.1 


9 = 0.5 


random SLHD (mean of 100) 


0.022 


0.027 


0.058 


random LHD (mean of 100) 


0.024 


0.026 


0.052 


Max. Entropy (0 = 0.05) SLHD 


0.016 


0.017 


0.025 


Max. Entropy(0 = 0.5) SLHD 


0.017 


0.018 


0.026 


Max. Entropy (9 = 1) SLHD 


0.019 


0.020 


0.027 


Maximin Distance SLHD 


0.020 


0.020 


0.028 



Table 10: Square root of MSE for Maximum Entropy SLHD, Maximin dis- 
tance SLHD, random LHD and SLHD over 400 randomly selected reference 
sites 

complicated and crucial to prediction accuracy. Sacks et al. (1989) suggested 
using maximum likelihood estimate of 9. However, the modeling process 
should not. be limited to Kriging. One advantage of using Latin hypercuhe 
designs is that they can facilitate almost any kind of model, parametric 
and non-parametric. Authors of this paper have used MARS (multivariate 
adaptive spline regression), GAM (generalized additive models) and second 
order polynomials to analyze computer experiments. 

The cyclone study was originally a case study in robust design. We choose 
this study to demonstrate the link between computer experiments and ro- 
bust designs. Robust design studies can also be carried out using computer 
models as presented by Welch et al. (1990). Orthogonal Latin Hypercube 
design and Symmetric Latin Hypercube design can be used in a robust, design 
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study as well. One can follow the response model approach of robust designs 
as proposed in Welch et al. (1990) and Shoemaker, Tsui and Wu (1991). 
First, establish a prediction model for both control and noise factors. Then, 
given the distribution of noise variables, estimate the variation of Y for each 
combination of control variables using the model obtained at the first stage. 
If a computer experiment is not expensive, one can skip the first step and es- 
timate the variation caused by the noise variable directly using the computer 
model. In that case, Latin hypercubes can serve as a sampling mechanism 
to obtain samples from noise variable, as it was first proposed by McKay et 
al (1979). 

5 Summary and Discussion 

This paper proposes a class of symmetric Latin hypercuhe designs (SLHDs), 
referred previously by Morris and Mitchell (1995) as "foldover designs", and a 
new columnwise-pairwise (CP) algorithm for searching optimal design within 
the SLHD class as well as within the regular LHD. 
We summarize the properties of SLHDs as follows. 

1. Thev are a good subset of LHDs with respect, to both entropy and 



maximin distance criteria (see Tables 2-4). 

2. As a generalization of Orthogonal Latin hypercube designs, SLHDs 
retain some orthogonality. The estimation of quadratic effects and bi- 
linear interactions is uncorrelated with the estimation of linear effects. 

3. The searching time of the CP algorithm is greatly reduced by restricting 
the search within the SLHDs (see Tables 5-8). The restriction does 
not significantly reduce the searching time of the simulated annealing 
algorithm, but it often leads to better designs (See Table 8). 

4. The global optimal LHD is not always a SLHD. Morris and Mitchell 
(1995 ) did an exhaustive search to find the optimal LHDs of small sizes. 
Not all of the true optimal designs they found are symmetric. 

Despite the fact that the true optimal LHDs do not. necessarily fall into the 
symmetric class. We recommend using the SLHD in computer experiments 
for two reasons. First, users will benefit from the orthogonal properties of 
SLHD as summarized above when they try to fit the data with a polynomial 
model. A non-symmetric LHD does not have such orthogonality. Second, 
as shown in Tables G and 8, by restricting the search within SLHDs, one 
could obtain approximately optimal designs in a more efficient, manner for 
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moderate to large-size designs. In fact, in these cases, an exhaustive search is 
usually prohibitive and one should be less concerned about whether a search 
method has the potential to reach the global optimum and more about how 
it can obtain a good design with reasonable computing effort. Especially for 
computer experiments, extra computing power could be spent on additional 
runs rather than obtaining a slightly better design. 

The performances of the three algorithms for searching optimal LHDs are 
summarized as follows: 

L The CP algorithm consistently outperforms the algorithm of Park (1994). 

2. For smaller designs, the CP algorithm is more efficient than the simu- 
lated annealing algorithm of Morris and Mitchell (1995). However, the 
latter can generate better large designs. 

One of the referees suggested that we briefly comment on the performance 
of optimal LHDs compare to other types of designs proposed for computer 
experiments in recent literature, such as orthogonal arrays, OA based Latin 
hypercubes and quasi-Monte Carlo lattices. The comparison of different kind 
of designs is one of the most important problems and deserves a thorough 
investigation that is beyond the scope of this paper. However, we would 



be glad to share some of our opinions. Unlike traditional designs for which 
the models are in known forms, the computer experimenter often has little 
idea which model in his/her statistical toolbox will best describe his/her 
complex computer model before an experiment is done and several kind of 
models are tried. Most of the proposed designs for computer experiments 
allow users to try many different models, linear or nonlinear, parametric or 
non-parametric. Among those, orthogonal arrays may not be appropriate for 
computer experiments since they do not take the advantage of flexibility of 
computer experiment in terms of changing levels. Their projections to low 
dimensions are only a few points so that the}' are not good for non-parametric 
regression methods. However, they are good for fitting low-order polynomial 
models. 

An optimal SLHD actually takes three criteria into consideration: the 
discrepancy of one-dimension projection optimized by the Latin hypercube 
structure, desired orthogonality inherited from the symmetric structure, and 
a third criterion (entropy or minimum distance) optimized through an algo- 
rithmic search. Therefore, we expect that optimal SLHD should perform very 
well with many modeling methods. Quasi-Monte Carlo lattice designs (Fang 
and Wang. 1994) are generated by some sequence which are asymptotically 



optimal in discrepancy measure. Since it spreads the design point evenly in 
the design space, it should have robust performance with different modeling 
methods. In particular, a glp (good lattice point) set is a Latin Hypercube. 
Bates et. al (1996) compared Latin hypercubes designs with lattice designs 
and found the quasi-Monte Carlo lattice design performed surprisingly well. 
Tang(1993) and Owen(1993) proposed a special type of Latin Hypercubes 
which are constructed based on orthogonal arrays. Such Latin hypercubes 
spread points evenly on ^-dimensional projections. The actual dimension 
t depends on the strength of the original orthogonal array. However, this 
approach only provides LHD at the sizes of which orthogonal arrays exist. 
In Table 11, we compare three LHDs. The first one is the fourth optimal 
SLHD listed in Table 6. The second one is a glp set of generating vector 
(25; 11 ,29,0,13). The third one is a LHD constructed based on 0.4(25; 5 4 ~ 2 ) 
using the procedure described in Tang (1993). It can be seen that in terms 
of entropy and minimum intersite distance, the optimal SLH is better than 
the glp and OA- based LH. The glp. however, is surprisingly good given the 
fact that it is easy to generate. Therefore, it could be a good choice if quick 
solutions to design problems in computer experiments are needed. On the 
other hand, the OA- based LH is far inferior to the other two designs. Since 
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Optimal SLH 



OA-based LH 



Minimum L\ 
Minimum Lo 



0.92(G) 
0.46(2) 
33.03 
21.68 
19.82 
15.96 
13.50 



0.75(24) 
0.41(12) 
39.09 
26.99 
24.97 
20.72 
17.9C 



0.54(13) 
0.29(13) 
50.17 
37.00 
34.77 
30.02 
26.89 



Entropy 0 = 0.05 



Entropj- 9 = 1 
Entropy 0 = 2 
Entropy 0 = 5 
Entropy 0 = 10 



Table 11: Comparison of three types of LHDs. 



a class of LHD can be constructed based on an orthogonal ana}', a similar 
algorithmic approach should be developed to find a better design within the 
class. 

We would like to see more research effort to compare those designs with 
respect to different perfor nicince measure. We think that a good design will 
not, necessarily score the highest for any particular criterion but will be rea- 
sonably high for all the criteria. 
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Abstract 

Traditional system dynamics studies rely heavily upon heuristics and experience. Never- 
theless, mathematical exploration techniques have been introduced as important elements 
for a successful study. We argue that the role of optimization in system dynamics studies is 
not to replace experience-based knowledge, but instead to augment, facilitate, and expand 
the heuristic exploration of a model. Accordingly, our approach involves narrowing the 
design space (using response surfaces) and the subsequent direct investigation of the simu- 
lation model (using heuristics). Response surfaces have received considerable attention in 
optimization because of their capability to replace complex models with analytic equations, 
thereby increasing computational efficiency. However, doubts exist as to the usefulness of a 
response-surface approximation of an approximation of reality (i.e., a system dynamics 
model). We demonstrate the usefulness of response surfaces in system dynamics studies with 
a case study involving a high-level model of an industrial ecosystem; our intent in using 
response surfaces is not to replace the simulation models with analytic equations, but instead 
to direct attention to regions within the design space of the original simulation with the 
most desirable performance. Recommended changes to a system are based directly on the 
simulation model, not on response surfaces, avoiding the added level of approximation 
inherent in response surfaces. The primary focus of the article is on the concept exploration 
approach, which is presented first. The case study towards the end is offered as supporting 
evidence. Copyright © 2000 John Wiley & Sons, Ltd. 

Sy$t. Dyn, Rev. 16, 75-90, (2000) 



Motivation 

System dynamics as a philosophy and feedback control modeling as an 
implementation of this philosophy are powerful in helping humans understand 
complex systems. The; approaches and methods associated with system dynam- 
ics have traditionally been heavily reliant upon the heuristic exploration of 
possible changes to a system. Nevertheless, optimization techniques have been 
introduced and used in system dynamics studies dating back to the early 1980s. 
Whereas some practitioners may have seen a future in which heuristics are 
completely abandoned in system dynamics studies, most of those exploring 
the role of optimization emphasized that mathematically rigorous optimization 
techniques should be coupled with heuristics to form a powerful approach 
for model improvement. Along these lines, our approach is a combination of 
mathematical tools and heuristics. 
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In particular, the mathematical tools we use are design of experiments, 
response surfaces, and a multi objective decision model. In using response sur- 
faces, a simulation model is approximated as an analytic equation in which 
the relevant responses are functions of significant control variables. Response 
surface techniques do not account for the fact that system dynamics models 
typically are rough approximations of reality — more so than a model that is 
analytic or based on extensive historical data. While system dynamics models 
can be constructed of systems that are analytic or based on extensive historical 
data, a strength of system dynamics lies in the ability to approximate the 
behavior and structure of poorly understood systems. The usefulness of a 
response surface approximation of an already highly approximated system 
dynamics model is questionable. Hence, our approach applies response sur- 
faces only to guide the search for an improved system; the final recommen- 
dations, however, are based directly on the simulation model. 

Additionally, our approach to 'optimization' does not involve finding one 
'optimal' solution, but instead finding a range of satisficing solutions. A sat- 
isfying solution is one that is 'good enough/ not necessarily optimal (Simon 
1996). A major reason for searching for satisficing solutions is that the optimal 
conditions in a model are most likely not optimal conditions in the real world, 
especially when modeling poorly understood systems. Therefore, by finding a 
range of solutions, more play between the model and reality can be accounted 
for while still achieving the desired responses. We begin by introducing pre- 
vious work on the use of optimization in system dynamics studies. 



Optimization and system dynamics modeling 

The inclusion of mathematical tools in the process of improving the behavior 
of system dynamics models has been examined by several people from varying 
perspectives. Multiple authors have worked in the area of optimizing the par- 
ameters of a model. An excellent starting point to learn about applying opti- 
mization to system dynamics models is clearly presented by Wolstenholm 
(1990) and Wolstenholm and Al-Alusi (1987). These authors do refer to their 
approach as 'heuristic optimization,' but it is unclear how heuristics influence 
their optimization process beyond the selection of objective functions and 
parameters to adjust. Keloharju and Wolstenholm (1989) emphasize the 
increased efficiency with which a parameter space can be explored using opti- 
mization techniques when compared to purely heuristic techniques. They do, 
however, maintain that human skill remains important in the formulation of 
the model. Interestingly, they find that several control variable combinations 



Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 



R. Bailey et al: Using Response Surfaces 77 



lead to very similar responses, an idea that supports the search for ranges of 
solutions introduced in our paper. Coyle's (1985) work is noteworthy in that 
he includes structural changes to a model through carefully formulating these 
structural changes with parameters. He, too, emphasizes the role of opti- 
mization as a tool complementary to the traditional heuristic approaches of 
system dynamic studies. 

Kleijnen (1995) includes design of experiments and response-surface meth- 
odology in his approach to optimizing the parameters of a model. Kleijnen intro- 
duces these tools as being 'objective, efficient, and effective* and applies them 
in a moderately standard manner. As in our paper, Kleijnen mixes heuristics 
and optimization to locate parameter values that produce good behaviors in the 
model The difference between our approach and Kleijnen's is that he narrows 
the design space with heuristics and applies response surfaces to optimize within 
this smaller region, whereas we apply response surfaces and optimization to 
narrow the design space and heuristics to locate a region of good behavior near 
the optimized solution. In addition, Kleijnen uses response surfaces to locate the 
final solution (and consequently checks the solution with the simulation model), 
while we advocate completely avoiding the added layer of approximation 
embedded in response surfaces when determining the final solution. 

For those not familiar with regression and design of experiments, Rotmans 
and Vrieze (1990) provide a straightforward review of these techniques, apply- 
ing them to a simulation model. Rotmans and Vrieze prefer the use of a central 
composite design, or CCD (which is also used in this paper), when applying 
design of experiments techniques to their simulation model. As an additional 
source, the use of orthogonal arrays versus Latin hypercubes in the sensitivity 
analysis of system dynamics models is discussed by Clemson et ah (1995). 
They point to response surfaces as being more appropriate when 'all parameter 
interactions need to be studied', which is often the case in optimization. 

One final class of approaches for enhancing the analysis and improvement 
of a system dynamics model is that of modal control, in which eigenvalues of 
the motion equations are used to synthesize new policy options (Macedo 1989). 
The main strength of using modal control theory is that new policy structures 
can be generated mathematically, as opposed to the heuristic methods of Coyle, 
Kleijnen, and Keloharju. Drawbacks of modal control theory include the necess- 
ary linearization of the model, the amount of computation, and the design 
of realistic policies from the synthetically generated policies (Mohapatra and 
Sharma 1985), Macedo has introduced a mixed approach in which modal con- 
trol and traditional optimization are sequentially applied in the improvement 
of a model. Though we do not focus on Macedo's work, his method appears to 
be promising. 
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The approach taken in this paper is similar to that of Keloharju in that the 
optimization of parameters is the primary focus. However, the implementation 
of design of experiments and response surface techniques draws upon the 
work of Kleijnen, A hidden assumption of Kleijnen, that a response-surface 
approximation of a simulation model (which is an approximation of reality) 
can be used to determine realistic options for parameter settings, is addressed 
and handled in this paper. Additionally, our focus is not on determining one 
'optimal* solution set of parameter settings, but instead on determining a range 
of solutions in which system behavior is acceptable. It is these two points of 
contrast to earlier research that form the contributions of this work. 

Frame of reference: robust concept exploration 

In recent years, design of experiments and response surfaces have gained in 
popularity as tools for aiding the exploration of concepts (Box and Draper 1987; 
Khuri and Cornell 1987; Myers and Montgomery 1995). The basic notion is to 
relate responses to controllable variables through analytic equations that are 
constructed through statistical methods. These analytic equations allow for 
quicker exploration of multiple concepts than do the simulation codes on 
which they are built. 

The robust concept exploration method (RCEM) is a formalization of the 
steps followed in using response surfaces in the design of a system (Chen et 
al, 1996). The RCEM is composed of four steps, starting with the identification 
of design parameters (Step l). Next, screening experiments are conducted, using 
design of experiments techniques to reduce the problem size (Step 2). Fol- 
lowing the screening experiments, response surfaces are created to map sig- 
nificant design parameters to response behaviors (Step 3). Finally, a 
multiobjective decision support problem (the compromise Decision Support 
Problem, or DSP) is formulated using the response surfaces and solved to facili- 
tate concept exploration and the identification of top-level design specifications 
(Step 4). 

The flow diagram in Figure 1 shows the flow of information during the RCEM 
process. Important to note is that final solutions are based directly on the 
response surface approximations of the original simulation code. 

When the simulation code is based on concrete relationships (such as those 
from physics or thermodynamics) or on extensive historical data, using 
response-surface approximations to generate final solutions has been shown to 
be effective (Koch et al 1996; Myers et al 1989; Simpson et al 1997). From a 
system dynamics perspective, however, models commonly are rougher 
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Fig. 1. Robust concept 
exploration method 
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approximations of reality — not based on concrete relationships, but instead 
on human interactions, aggregated trends, and less precise information. The 
usefulness of a response-surface approximation of a highly approximated simu- 
lation model has not been directly addressed in the literature. Our approach to 
using response surfaces with system dynamics emerges from this dilemma. 



Our approach 

A response-surface approximation of a simulation model that is already an 
approximation of reality leads to two layers of approximation between the 
solutions and the real world. If possible, there is a need to eliminate one level 
of approximation. A way to remove one level of approximation is by adding a 
loop to the baseline RCEM. With this approach, shown in Figure 2, response 
surfaces and a compromise DSP are used as attention-directing tools, forming 
an information feedback loop (Bailey 1997). 
Instead of taking the solutions directly from a compromise DSP, the results 
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from a compromise DSP (which are based on response-surface approximations 
of the simulation model) are fed back to the simulation model. The regions 
near the compromise DSP solutions are explored with the simulation model in 
search of a range of acceptable solutions. The response surfaces are still 
extremely useful in identifying a promising region of the design space, but not 
in identifying final solutions. The process of creating the response surfaces is 
termed the comprehension loop, because the relationships between responses 
and control variables are identified and learned in this loop. 

The added step in the augmented RCEM (shown in Figure 2) involves return- 
ing to the simulation program and using the solutions from the compromise 
DSP formulations as guides to explore the design space. Ranges of control 
variables for further exploration are determined with two heuristic principles 
(Bailey 1997): 

1. If a compromise DSP solution for a control variable is not on a bound, then 
two points — one to either side of the compromise DSP solution — are selected 
for further exploration. 

2. If a compromise DSP solution for a control variable is on a bound, then 
instead of exploring one point on each size of the compromise DSP solution, 
two points are explored on the side of the compromise DSP solution within 
the bound. 

The points used for further analysis (called exploration points), therefore, sur- 
round the compromise DSP solution in orthogonal directions. Values for the 
deviation variables (i.e., the deviations from the desired responses) are found at 
the exploration points with the original simulation program (not the response- 
surface models) and are then compared to the compromise DSP, or baseline, 
solution. Consider the example shown in Figure 3, in which, for simplicity, a 
two-dimensional design space is explored (in the case study, the design space 
is ten-dimensional). 

The compromise DSP solution is found and input into Step 1. Points orthog- 
onal to the compromise DSP solution are found for each of the two variables 
in Step 2. In Step 3, the deviation variables at each one of the solution explo- 
ration points are determined with the simulation model (these deviations are 
shown in Figure 3 as vertical columns). Finally, in Step 4, these deviation 
variables are compared to the deviation variables for the compromise DSP 
solution. If the deviation variables are within a certain percentage of the 
compromise DSP solution, then that point is included in the solution range. A 
designer must determine the allowable amount of fluctuation for each deviation 
variable (in our case study, a ten percent fluctuation is allowed). If the responses 
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Fig. 3. Solution 
exploration process 
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are for a very precise system (e.g., a laser used in surgery), a smaller percentage 
fluctuation is appropriate. 
We discuss two main points concerning the solution exploration strategy: 

1. The exploration points are determined heuristically. 

2. The solution ranges are based on main effects (Bailey, 1997). 

The heuristic determination of the exploration points leads to more flexibility 
in the solution exploration scheme. A designer has considerable knowledge 
concerning the sensitivity of the control variables after formulating a decision 
support problem, developing a simulation model, building response surfaces, 
and analyzing results from compromise DSPs. Setting concrete rules for deter- 
mination of the exploration points would limit the process more than enhance 
it. 

The sole use of main effects (i.e., first-order effects) in developing ranged sets 
of solutions is based on the premise that the interactions and second-order 
effects do not have much influence on the behavior of the system when small 
perturbations are made in the values of the control variables. Large changes in 
the values of the control variables are not made in the solution exploration 
process, because interaction effects can begin to influence the system behavior 
more strongly. Quadratic effects, however, are modeled by the response sur- 
faces earlier in the concept exploration process. This is necessary to obtain 
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response surfaces that approximate the behavior of the simulation model accu- 
rately enough to search the entire design space. 

A benefit of exploring the regions surrounding each control variable is the 
development of a ranged set of solutions, as opposed to a singular, point solu- 
tion. With the augmentation made to the RCEM, design freedom is maintained 
while design knowledge is increased with the solution exploration process. 
The use of response surfaces and the compromise DSP as attention-directing 
tools (instead of solution-determining tools) is a significant step to maintaining 
design freedom while gaining information about the system. Instead of using 
these tools to narrow the design space to one final solution, response surfaces 
and the compromise DSP are used to expand the information about the design 
space, which is narrowed to a final range of solutions after an additional step. 

The simulation slot in the augmented RCEM (see Figure 2) is where domain- 
dependent information enters the concept exploration process. The simulation 
model presented in this paper is of a system of industries in Kalundborg, 
Denmark. 

Case study: a system of industries in Kalundborg, 
Denmark 

Over the past 30 years in Kalundborg, Denmark, a system of industries has 
developed in which by-products and wastes from one company are used by 
other companies in the system. For instance, industrial gypsum collects in the 
scrubbers at Asnacs, a coal and gas burning power plant. This industrial gypsum 
is sold to Gyproc, where it is used instead of mined gypsum to make wallboard. 
A series of such interactions exists in Kalundborg and is described in more 
detail in Bailey et ah (1997), Ehrenfeld and Gertler (1997) and Frenay (1995). 

A simulation model of a section of Kalundborg (including the industries with 
the most interactions with other industries) was built, validated, and used to 
generate data in STELLA (STELLA is described in Chichakly et ai 1994), While 
the model and its validation are described in detail by Bailey (1997), only a 
brief description is presented here. In Figure 4, a simplified version of the flow 
diagram of the simulation model is shown. 

There are three primary accumulations (steam, gas, and wallboard) that are 
each individually part of balancing loops. The targets for gas and wallboard 
are exogenous variables, while the target for the steam accumulation is set by 
current usage of steam (steam cannot be inventoried indefinitely — it is per- 
ishable — so current accumulation levels must be based on steam usage). The 
interactions among the three companies connect the three balancing loops. For 
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Fig. 4. Simplified flow 
diagram of Kalundborg 
model 
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instance, since steam is used in the purification of gas, an information arrow 
connects the gas inflow to the steam outflow. 

The purpose of the concept exploration process is to design a system that is 
robust to a disturbance in the system. The primary responses used to evaluate 
system behaviour are based on the three accumulations, while the controllable 
variables deal with how the accumulations respond to a disturbance. For the 
purposes of this paper, only one disturbance— namely, a step change in steam 
demand — is presented. The responses are selected to characterize the rate of 
recovery, mean value, and standard deviations of the accumulations as they 
attempt to reach the stepped target value. Additionally, the total amount of 
waste generated in the system is calculated and used as a response. Altogether, 
there are eight responses used to characterize the robustness of the system and 
the waste generated. 

The three types of control factors occur several times within the model, 
resulting in a total of ten control variables. These delay times, averaging times, 
and correction times each have physical significance in the actual system: 

• Delay time: time between the occurrence of a discrepancy between the goal 
and actual level of an accumulation and the occurrence of first reaction. 
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Fig. 5. (a) Face-centered 
and (b) inscribed 
central composite 
designs 
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♦ Averaging time: time over which flows are averaged; averages of production 
flows are used in determining present production rates. 

• Correction time: time over which a discrepancy between the goal and actual 
levels of an accumulation is intended to be corrected. 

All of the decisions in the simulation model relate to changing the production 
rates of steam, gas or wallboard to maintain the desired accumulation levels. 
Additionally, the decisions are parametric in nature; they do not involve struc- 
tural changes in the model. The ten control variables are related mathematically 
to the eight responses through response surfaces. The basic structure of the 
response surface formulation is presented in the next section. 

Design of experiments and response surface formulation 

To develop quadratic response surfaces, experiments must be conducted with 
the simulation model to determine the effects of each control variable and the 
interactions between control variables. Among the various types of exper- 
imental designs for fitting such response surface models, the central composite 
design (CCD) is probably the most widely used (Montgomery 1991). For each 
problem, a slightly different strategy will need to be used to fit response 
surfaces. The strategy we used is outlined here and is described in more detail 
by Bailey (1997). 

First, a face-centered CCD is used to relate the ten control variables to the 
eight responses and screen off insignificant factors. A three-factor face-centered 
CCD is shown in Figure 5(a) with experimental points at the corners (i.e., the 
bounds for the three factors), the centers of the faces, and in the center of the 
entire space. In these screening experiments, we found that only four of the 
eight responses are affected by the step change in steam demand disturbance. 
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Additionally, the number of significant control variables is reduced from ten 
to seven. 

With the problem reduced in size, a second experiment, using an inscribed 
CCD, is performed. A three-factor inscribed CCD, shown in Figure 5(b), uses 
the same number of points as a face-centered design. The difference is that the 
corner points are moved from the extreme bounds inward, while the face- 
centered points remain in the same positions. The response surfaces obtained 
from this second round of experiments are used in the compromise DSP to 
search for promising regions in the design space (as defined by the values of 
the control variables). The response surfaces are judged by the R 2 value (coef- 
ficient of determination), which represents the proportion of the total variation 
in the response than can be accounted for by the quadratic relationship with 
the control variables (Simpson et al. 1997; Walpole and Myers 1989), An R 2 of 
zero means that none of the variation in response can be accounted for by the 
response surface, whereas an R 2 of one means that all of the variation is 
accounted for by the response surface. In our case study, all response surfaces 
have an R 2 of 0.91 or higher. 

Compromise DSP solution 

The compromise DSP is the mathematical construct through which the various 
metrics and tools of the RCEM are integrated. It is a multiobjective decision 
model, which is a hybrid formulation based on mathematical programming and 
goal programming used to find a set of system variables that satisfy system 
constraints while achieving a set of conflicting goals as well as possible (Mistree 
et aL 1993). The system descriptors, namely system and deviation variables, 
system constraints, system goals, bounds, and the deviation function, are 
described in detail by Mistree et al. (1993). The mathematical formulation for 
a compromise DSP is as follows: 

Given: an alternative to be improved. Assumptions used to model the domain 
of interest. 
The system parameters: 

n number of system variables 

q inequality constraints 

p+q number of system constraints 

m number of system goals 

gi(X) system constraint function 

f k {d t ) function of deviation variables minimized at priority level J: for 
the preemptive case. 



Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 



86 System Dynamics Review Volume 16 Number 2 Summer 2000 



Find: control variables and deviation variables 
Xi i = 1, . . n;dr,dji= 1, ...,m. 

Satisfy: 

System constraints (linear, nonlinear) 

gi (X)=0;i=l,..,p 

g(X)^ 0;i = p+l,..,p+q 
System goals (linear, nonlinear) 

Aj(X )+dr-dt = G i ;i=l,...,m 
Bounds 

XT in ^X, -,^X? ax ;i=l n 

dr,d+ ^ 0;i=l,...,m 

dr-d7 = 0;i=l,...,m 

Minimize: 

deviation function 

z=[f,(dr,dt) fjdr.dw 

The response surfaces are used in formulating the goals in the compromise 
formualtion. The particularization of the compromise DSP for the step change 
in steam usage is presented by Bailey 1997; Bailey et ah 1998. 

Solution exploration 

Following the procedure outlined in the section above describing our approach, 
solution exploration points are identified that surround the compromise DSP 
solution. At each of these points, the four deviation variables are calculated 
using the simulation model of Kalundborg (not the response surfaces). These 
values of the deviation variables are then compared to the values of the devi- 
ation variables for the compromise DSP solution. In Table 1, to illustrate this 
process, the solution exploration points and associated deviation values are 
shown for three of the ten control variables. 

A solution exploration point is accepted into the final range of specifications 
only if none of the four values of the deviation variables increases to more than 
ten percent of the values of the deviation variables from the compromise DSP. 
In Table 1, all but two deviation variables for the solution exploration points 
are within an acceptable range. The two unacceptable points (the steam cor- 
rection time, STMCT, at 3.5 months and the steam delay, STMDEL, at 1.0 
months) are those associated with runs #2 and #4. All exploration points that 
have values of the deviation variables within ten percent of the compromise 
DSP values (runs 1, 3, 5 and 6) are included in the final solution ranges. 
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Table 1. Solution 

exploration step Control variables (months) Deviation variables 



Run STMCT STMDEL ... ATUSE d, d 2 d 3 d 4 



cDSP 


3.25 


0.5 






12 


0.250 


1. 


85e- 


■3 


0.0326 


0.945 


1 


3 


0.5 






12 


0.219 


1. 


85e- 


■3 


0.0337 


0.948 




Percentage 


difference 


from 


cDSP 


solution 


-12.4 


0 






+3.37 


+3.17 


2 


3.5 


0.5 






12 


0.283 


1. 


85e- 


■3 


0.0317 


0.942 




Percentage 


difference 


from 


cDSP 


solution 


+13.2 


0 






-2.76 


-0.317 


3 


3.25 


0.75 






12 


0.263 


1, 


85e- 


■3 


0.0341 


0.944 




Percentage 


difference 


from 


cDSP 


solution 


+5.2 


0 






+4.6 


-0.106 


4 


3.25 


1.0 






12 


0.276 


1. 


85e- 


-3 


0.0356 


0.944 




Percentage 


difference 


from 


cDSP 


solution 


+10.4 


0 






+9.20 


-0.106 


5 


3.25 


0.5 






10 


0.272 


1. 


90e- 


■3 


0.0351 


0.952 




Percentage 


difference 


from 


cDSP 


solution 


+8.8 


+2. 


,70 




+7,67 


+0.741 


6 


3.25 


0.5 






11 


0.261 


1. 


.87e- 


■3 


0.0338 


0.948 




Percentage 


difference 


from 


cDSP 


solution 


+4.4 


+1, 


.08 




+3.68 


+0.317 



* Control valiables not shown are set to cDSP solution variables for the runs shown 
STMCT: steam correction time; STMDEL: steam delay; ATUSE: steam usage averaging time 



Table 2. Ranged set of 

solutions (in months) Desi S n s P ace bounds 



Control variable Upper Lower Solution values 



Steam correction time 


0.5 


6.0 


3.0-3.25 


Steam delay time 


0.5 


2.0 


0.5-0.75 


Steam generation averaging time 


0.5 


12.0 


5.0-6.0 


Steam usage averaging time 


0.5 


12.0 


10.0-12.0 


Gas correction time 


0.5 


6.0 


0.5 


Gas delay time 


0.5 


2.0 


0.5-1.0 


Gas purification averaging time 


0.5 


12.0 


9,25-11.25 


Wallboard correction time 


0.5 


6.0 


4.0-6.0 


Wallboard delay time 


0.5 


2,0 


0.5-1.0 


Wallboard production averaging time 


0.5 


12.0 


10.0-12.0 



The range of solutions obtained for the case study, developed using the 
solution exploration process, is shown in Table 2. The solutions uniformly 
encourage the use of short delay times and long averaging times, while cor- 
rection times display no general trend. Shorter delay times are representative 
of situations in which decisions are made with more current information. Inter- 
esting to note is that the solution exploration process has provided a range of 
solutions that are satisficing (the delay times do not have to be as short as 
possible to achieve desirable responses). Long averaging times are characteristic 
of systems that react to disturbances gradually (not abruptly). Therefore, a long- 
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term approach to handling disturbances in the Kalundborg model results in 
better performance. 

The verification of these results is not shown here; however, it may be found 
in its entirety in Bailey (1997) and in an abbreviated form in Bailey et ah (1998). 
Through the verification process, the solution ranges are shown to produce the 
desired robust behaviors in the model. 



Closure 

There is certainly an opportunity for using design of experiments and response 
surfaces in the improvement process of a model. The transferal of these tools 
from other disciplines, however, must be done with care — assumptions con- 
cerning their use must be examined. In the case where a model is a rough 
approximation of reality, the use of response surfaces to determine the final 
solution is questionable since these response surfaces are further approxi- 
mations. Nevertheless, response surfaces are powerful and have potential for 
increasing the efficiency and effectiveness of system dynamics studies. Our 
approach uses response surfaces to direct attention to a particular location in 
the design space, and then — returning to the simulation model — develops a 
range of solutions that will lead to the desired behaviors. The particular 
implementation of this approach, based on the robust concept exploration 
method, is outlined in this paper and explored through an example. 

One potential drawback to this approach is that the added step — returning 
to the simulation model — will decrease the efficiency of the process. While 
some efficiency is lost, the total number of simulation runs added is only twice 
the number of control parameters — in our example, 20 runs. Including all of 
the data entry and analysis, the added step required roughly two hours to 
perform for the Kalundborg example. The benefits of the solution exploration 
step — using optimization to find a range of solutions that are based directly on 
the simulation model — outweigh this slight loss in efficiency. 

Approaches such as the one presented here, in which heuristics and math- 
ematical tools are combined in a system dynamics study, are developing as 
important tools in the improvement of model behavior. For an experienced 
person, the emphasis might lean towards heuristics, whereas for a novice the 
mathematical tools can help that person gain experience while reducing the 
number of major mistakes made along the way. Each component, the heuristic 
and the mathematical, is an important element to a successful study. 
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Sequential approach, Weibull. 

Summary & Conclusions — This tutorial explains statistically 
designed experiments which provide a proactive means to improve 
reliability as advocated by Genichi TaguchJ. That is, by systematic 
experimentation, the important parameters (factors) affecting 
reliability can be identified along with parameter values that yield 
reliability gains. In addition to improving reliability, Taguchi's 
robust design can be used to achieve robust reliability; that is, to 
make a process or product reliability insensitive to factors which 
are hard or Impossible to control. Robust design is also implemented 
using statistically designed experiments. This paper presents classes 
of experimental plans for reliability improvement and robust 
reliability. An important feature of the reliability data collected 
from such experiments is censoring which occurs when some of 
the experimental units have not failed by the end of the experi- 
ment. Consequently, the analysis methodology must account for 
these censored data which are likely to occur in light of the ever 
increasing reliability of today's products. Several appropriate 
methods are discussed briefly. These experimental plans and 
analysis methods are illustrated using three documented ex- 
periments which improved fluorescent lamp and industrial ther- 
mostat reliability and which achieved robust reliability for night- 
vision goggles. 



1. INTRODUCTION 

Statistically designed experiments have been used exten- 
sively for estimating or demonstrating existing reliability, eg, 
most of the examples in Nelson [15]. Until recently, designed 
experiments appear to have seldom been used to improve 
reliability by identifying the important parameters (factors) af- 
fecting reliability out of many potentially important ones. For 
example, Taguchi [25 - 27] advocated using designed ex- 
periments as a proactive means for improving reliability and 
provided examples of improving clutch springs and fluorescent 
lamps. Taguchi is perhaps best known for robust-design, whose 
aim is to make processes/products insensitive to noise factors 
which are hard or impossible to control. Such products/pro- 
cesses are robust to the noise factors. Examples of noise fac- 
tors include manufacturing variables that cannot easily be con- 
trolled, and environmental conditions in which the product is 
used. This important paradigm for improving products/pro- 



cesses, which attracted the attention of industry in the 1980s 
[11] can also be applied to reliability. To ensure good stability 
and adequate reliability, Taguchi [25: page 149] recommend- 
ed that noise factors be considered in any experiment to im- 
prove reliability when it is practical to do so. 

Since the early 1980s, several experiments for improving 
reliability have been documented. Specht [23] reported the im- 
provement of heat-exchanger reliability in a commercial heating 
system; Montmarquet [14] discussed the improvement of drill- 
bit reliability in a multilayer printed circuit board drilling opera- 
tion; and Reed [19] presented the improvement of night-vision 
goggle reliability. From an early application of Taguchi's 
methodology at AT&T, Phadke [16] discussed the improvement 
of router-bit reliability in a printed circuit board cutting opera- 
tion. Condra [4] gave several examples from the electronics in- 
dustry. Taguchi* s robust design philosophy figures prominent- 
ly in Condra's book. Recently, Bullington, Lovin, Miller, 
Woodall [3] reported on the improved reliability of industrial 
thermostats. 

This tutorial illustrates how statistically designed ex- 
periments can be used to improve reliability and to achieve 
robust reliability. Classes of experimental plans for improving 
reliability and achieving robust reliability are discussed with 
examples in sections 2 & 3, respectively. Section 4 briefly 
discusses analysis methodology for extracting the information 
from the experimental data. The methodology must account for 
censoring, an important feature of such experiments, which oc- 
curs when some of the experimental units have not failed by 
the end of the experiment; this type of censoring produces type-I 
right-censored data. Other types of censoring typically arise in 
such experiments. Where units cannot be monitored continuous- 
ly, units must be inspected periodically until failure. Periodic 
inspection produces left-censored data for units failing before 
the first inspection and interval -censored data otherwise. Ap- 
propriate analysis techniques which handle censored data are 
illustrated in sections 5-7 which analyze experiments to im- 
prove fluorescent-lamp and thermostat reliability and to achieve 
robust reliability of night-vision goggles. 

Acronyms 1 

ML maximum likelihood 
MLE ML estimate. 

Assumptions 

1. Lifetimes are j- independent and either Weibull or 
lognormally distributed. 

2. Some lifetimes are censored. 



'The singular & plural of an acronym are always spelled the same. 
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Notation 
2 k ~ p . 

3 k ~ p [2-level, 3-levelJ fractional factorial, k factors 
Yi response i 
Tj lifetime / 

Xj vector / of covariate values associated with Yi 
f} vector of regression parameters or effects 
a scale parameter 
€ t error i associated with Y ( 

Control vector of control -factor levels, also covariate values 
in the regression model involving only the control 
factors 

*noise vector of noise-factor levels, also covariate values in 
the regression model involving only the noise factors 

^control x noise covariate values in the regression model cor- 
responding to control-noise interactions 

{i€CENS}, {/€FAIL} set of [censored, failure] data 

<£(•) likelihood function 

Aj level j of factor A. 

Other standard notation is given in ' 'Information for Readers 
& Authors* 1 at the rear of each issue. 



2. EXPERIMENTS FOR IMPROVING RELIABILITY 

While there are potentially many factors (parameters) 
that affect reliability, some factors are more important, viz, 
have a bigger impact on reliability as the values of these 
factors are changed. These important factors can be identified 
empirically through experimentation which involves making 
deliberate changes in the factor values (levels) and observing 
the resulting reliability. Besides identifying the important fac- 
tors, levels for these factors that yield reliability gains can 
be recommended. Statistically designed experiments provide 
a systematic & efficient plan of experimentation to achieve 
these goals so that several factors can be studied simultaneous- 
ly. Designed experiments, which have successfully been used 
to improve other quality characteristics can also be used to 
improve reliability. 

2.1 Full Factorial Designs 

Some terminology is helpful in describing the plans. The 
plan of experimentation is the experimental design, or simply 
design. The experimental design consists of a list of runs 
(a run is a combination of levels at which the factors in the 
experiment are set). The number of runs in the experimental 
design is the run-size. The experiment then involves making 
units according to the conditions specified by the runs in the 
experimental design and lifetesting these units to failure. Table 
1 gives an experimental design for 3 factors (A, B, C) with 
each factor being studied at 2 levels (1, 2). Run 1 indicates 
that units are made with all the factors set at their respective 
level Vs. This particular design is a '2-level full-factorial' 
since it consists of all possible combinations of the 2 levels 
for the 3 factors. 



TABLE 1 

2-Level Full Factorial Design for 3 Factors 



Factor 



Ran ABC 



1111 

2 112 

3 12 1 

4 12 2 

5 2 11 

6 2 12 

7 2 2 1 

8 2 2 2 



The run size for a 2 -level full -factorial design in k factors 
is 2*, which quickly becomes prohibitive for more than 5 
factors. Designs with more than 2 levels allow curvature, 
interactions, or nonlinear effects to be assessed but require 
more runs. Even for 3 factors, the run size of the full-factorial 
design is 3 3 =27. Because of their large run sizes, full- 
factorial designs tend not to be used in an initial experiment 
unless there are only a few potentially important factors to 
be studied. 

2.2 2-Level Fractional Factorial Design 

For the typical industrial situation, many factors need 
to be studied in a relatively few runs. A sequential approach 
to experimentation provides one such strategy. An initial ex- 
periment using only a few levels (often 2) for each factor 
is performed to screen out the unimportant factors. A follow- 
up experiment involving much fewer factors but at more levels 
can then be used to explore the response-factor relationship 
in more detail. For the initial experiment, a subset or fraction 
of the full factorial design (a fractional factorial) can be used. 
There are two types of 2-level designs: 

• geometric (2*~ p ) designs [2] 

• non-geometric designs [17]. 

The notation for the geometric designs indicate the degree 
of fractionation, viz* a 2~ p fraction of a full factorial with 
run size 2 k ~ p . 

Taguchi [26: page 930] provided an example of a 2 k ~ p 
design which was used to improve the reliability of fluores- 
cent lamps. The experiment studied five 2-level factors (A 

- E) in 8 runs using a 2 5-2 design (a quarter of a full fac- 
torial) as given in table 2. Eight types of lamps were built 

— contrasted with a full-factorial design which would have 
required 32 types. No further details on the factor names 
and levels were provided. Two lamps were made at each 
run, and life testing was conducted over 20 days with inspec- 
tions for failure every 2 days. Seven of the 16 lamps had 
not failed by the 20-day inspection, which yielded right- 
censored data. 

Bullington ex al [3] provided an example of a 12-run 
Plackett-Burman design which was used for an experiment 
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TABLE 2 
Design & Lifetime Data 
[Fluorescent-Lamp Experiment] 



Factor 



A 


B 


C 


D 


E 


? 


7 


Lifetime* 


1 


1 


1 


1 


1 


1 


1 


(14,16) 


(20, oo) 


1 


I 


2 


2 


2 


1 


2 


(18,20) 


(20,oo) 


1 


2 


1 


1 


2 


2 


2 


(08.10) 


(10,12) 


1 


2 


2 


2 




2 


1 


(18,20) 


(20,oo) 


2 


1 


I 


2 


1 


2 


1 


(20,oo) 


(20,oo) 


2 


1 


2 


1 


2 


2 


2 


(12,14) 


(20,oo) 


2 


2 


1 


2 


2 


I 


2 


(16,18) 


(20,oo) 


2 


2 


2 


1 


I 


1 


1 


(12,14) 


(14,16) 



•The lifetime is uncertain. The 2 endpoints of the interval, within which the 

lifetime lies, are shown in 0. An end point of "oo" implies that the item had 

not failed at the end of the test. The time-unit is days. 

?-This design could have accomodated 2 more factors in these 2 unused 

columns. 



to improve the reliability of industrial thermostats. Eleven fac- 
tors (A - K) were studied using the design in table 3 in which 
10 thermostats were manufactured at each of the 12 run set- 
tings, viz, 12 types of thermostats were built. (A full-factorial 
design would have required 2 11 = 2048 types.) These factors 
were chosen from many across a 14-stage manufacturing pro- 
cess and include: E — beryllium-copper grain size, H — heat 
treatment, J — power element electroclean, and K — power 
element plating rinse. Each factor was studied at 2 levels, eg, 
E — grain size of 0.008 & 0.018 inches, and H — 45 & 240 
minutes. Table 3 also presents the lifetime data; the experiment 
was stopped at 7432 k-cycles resulting in 22 right-censored 
observations at runs 1, 6, 11. 

While highly fractionated 2 k ~ p and Plackett-Burman 
designs are ideally used as screening designs, in practice, the 



initial experiment might be the only one performed. Consequent- 
ly, a properly chosen 2 k ~ p design can allow some interactions 
to be studied. For example in the fluorescent-lamp experiment, 
besides the factors A - E main effects, the experimenter also 
thought that the A x B interaction might be important; this ex- 
perimental design (table 2) allows the Ax B interaction to be 
estimated. Factor main-effects refers to the additive effects of 
the factors on reliability. The interaction between two factors 
indicates the degree of non-additivity of the factor effects; ie t 
if an interaction is present, the effect on reliability of changing 
the levels of one factor depends on the level of the other factor. 
Consequently, the presence of an interaction can impact the 
recornmendations of levels for setting the important factors. See 
[2] for more discussion. For Plackett-Burman designs, [9] show- 
ed that some information on interactions can be obtained. Ref 
[9] proposed a sequential strategy which first identifies the im- 
portant factor main effects. Then, interactions involving these 
identified factors are entertained for possible importance. 

2.3 Multi-Level Designs 

Taguchi [26] often initially uses designs with more than 
2 levels. These include the 3 k ~ p designs and mixed-level 
designs such as the 18-run design which can be used to study 
one 2-level factor and up to seven 3-level factors. Refs [5, 
28] described other mixed-level designs. An example of the 
use of these designs is the clutch-spring experiment in Taguchi 
[25: chapter 9] which used a 3 7 " 4 design to study 7 factors 
in 27 runs, a 1/81 fraction of a 3-level full factorial design. 
Thus, each of the factors was studied at 3 levels. Phadke 
[16] also used a mixed-level 32-run design to study two 4-level 
factors and seven 2-level factors in a routing process. One 
of the 4-level factors was router bit type so that 4 types of 
bits were studied. Thus, qualitative factors such as the router 
bit type can be investigated as well as quantitative factors 
such as temperature. 



TABLE 3 
Design & Lifetime Data 
[Thermostat Experiment] 



Design 



A 


B 


C 


D 


E 


F 


G 


H 


I 


J 


K 






Ordered Lifetime Data (k-cycles) 










1 


1 


1 


1 


1 


1 


1 


1 


1 


1 


957 


2846 




















1 


1 


1 


1 


2 


2 


2 


2 


2 


2 


206 


284 


296 


305 


313 


343 


364 


420 


422 


543 




1 


2 


2 


2 


] 


1 


1 


2 


2 


2 


63 


113 


129 


138 


149 


153 


217 


272 


311 


402 




2 


1 


2 


2 


1 


2 


2 


1 


1 


2 


76 


104 


113 


234 


270 


364 


398 


481 


517 


611 




2 


2 


1 


2 


2 


1 


2 


1 


2 


1 


92 


126 


245 


250 


390 


390 


479 


487 


533 


573 




2 


2 


2 


1 


2 


2 


1 


2 




1 


490 


971 


1615 


6768 














2 


1 


2 


2 


1 




2 


2 




2 


1 


232 


326 


326 


351 


372 


446 


459 


590 


597 


732 


2 


1 


2 


1 


2 


2 


2 


1 


1 


1 


2 


56 


71 


92 


104 


126 


156 


161 


167 


216 


263 


2 


1 


1 


2 


2 


2 


1 


2 


2 




1 


142 


142 


238 


247 


310 


318 


420 


482 


663 


672 


2 


2 


2 


1 


] 




1 


2 


2 


1 


2 


259 


266 


306 


337 


347 


368 


372 


426 


451 


510 


2 


2 


1 


2 


1 


2 


1 


1 


1 


2 


2 


381 


420 


















2 


2 


1 


I 


2 




2 


1 


2 


2 


1 


56 


62 


92 


104 


113 


121 


164 


232 


258 


731 



•Test was stopped at 7342 (right censored); this and subsequent items did not fail. 



HAMADA: SDE TO IMPROVE RELIABILITY AND ACHIEVE ROBUST RELIABILITY 



209 



In contrast with the initial use of multi-level factor designs, 
the sequential approach to experimentation uses multi-level 
designs in a follow-up experiment. Box & Draper [1] gave 
response-surface designs which allow the response-factor rela- 
tionship to be explored in more detail. For example, an inter- 
nal document from a North American automobile manufac- 
turer reported the use of an 8-run experiment to screen 7 
factors. (The same design given in table 2 consisting of the 
first 7 columns was used.) Four factors were identified, and 
were studied further using the 27-run Box-Behnken design 
in table 4. (-1, 0, +1) denotes the 3 levels for each factor; 
each of the groups of rows, 1-2, 4-5, 7-8, specify 
4 runs since the ±1 notation means that all combinations 
of levels 1 & 3 are used. Other response surface designs 
such as the central-composite designs (with 5 factor-levels) 
can be employed. See [1] for more details. 

TABLE 4 

Box-Behnken Response-Surface Design for 4 Factors 



Factor 



A 


B 


c 


D 


±1 


±1 


0 


0 


0 


0 


±1 


±1 


0 


0 


0 


0 


±1 


0 


0 


±1 


0 


±1 


±1 


0 


0 


0 


0 


0 


±1 


0 


±1 


0 


0 


±1 


0 


±1 


0 


0 


0 


0 



Nomenclature 

INT an intercept 

A,B,C factor main effects 

A x B t A x C,B x C 2-factor interactions. 

To obtain the main effects values, [1 , 2] in table 1 are replaced 
by [-1, +1]. Values for the interactions are obtained by 
multiplying those from the corresponding main effects, eg y the 
A x# values are obtained by multiplying the A & B main-effect 
values. Ref [13J discussed covariate values for more than 2 
level-factors. 

TABLE 5 

Covariates for 2-Level 3-Factor Full-Factorial Design 



Values 



Run 


INT 


A 


a 


c 


AxB 


AXC 


BxC 


1 


1 


-1 


-i 


-l 


+ 1 


+ 1 


+ 1 


2 


1 


-1 


-t 


+ i 


+ 1 


-1 


-1 


3 


1 


-1 


+i 


-l 


-1 


+ 1 


-I 


4 


1 


-1 


+i 


+ i 


-1 


-1 


+ 1 


5 


1 


+ 1 


-l 


-l 


-1 


-1 


+ 1 


6 


I 


+ 1 


-l 


+ i 


-1 


+ 1 


-1 


7 


I 


+ 1 


+ i 


-i 


+ 1 


-1 




8 


1 


+ 1 


+ i 


+ i 


+ 1 


+ 1 


+ 1 



Section 4 discusses analysis methodology for fitting these 
models in (1) using censored data. Sections 5 & 6 present 
analyses of the fluorescent-lamp- and thermostat-experiments 
respectively. 



2.4 Analysis of Lifetime Data 

The lifetime data from these experimental designs can be 
analyzed using a parametric regression model such as the lognor- 
mal or Weibull. These models have the form [12]: 

Yt = log(7;-) = x?-p + a-€ if i'=l,..., n. (1) 

For the Weibull model, the errors {ej are i.i.d. standard 
extreme- value r. v, , whose Sf is exp(exp(-u>)). For the lognor- 
mal model, the errors {e,} are i.i.d. standard s-normal r.v., 
whose Sf is gaufc(w). These models assume that the factors af- 
fect the mean Iog(lifetime) (through x T -fi) and do not affect the 
common scale-parameter a. In terms of these models, the 
analysis first seeks to identify the non-zero components of the 
vector of effects 0. The recommended levels for setting the im- 
portant factors are those levels which maximize the mean 
log(Iifetime) x T -$. 

2.5 Design vs Corresponding Covariates 

Table 5 gives the corresponding covariate values for the 
2-level 3-factor full-factorial from table 1. 



3. EXPERIMENTS FOR ACHIEVING ROBUST 
RELIABILITY 

Taguchi's robust design is also referred to as parameter- 
design because its objective is to find levels of engineering 
parameters (control factors) that yield a robust product/process, 
ie % that make the product/process insensitive to the variation 
of hard- or impossible-to-control noise factors. Robust design 
is therefore different from the approach of handling sources of 
variation by control, which can be cosdy. For example, the Ina 
Tile Company was faced with reducing an unacceptable amount 
of variation in their tile-size caused by an uneven temperature 
distribution in the kiln [11]. Rather than purchasing an expen- 
sive kiln which would have controlled the temperature distribu- 
tion better, it was found through designed experiments that in- 
creasing the lime content in the tile formulation decreased the 
tile size variation by a factor of 10. In other words, a tile for- 
mulation was found that was insensitive to the existing oven's 
uneven temperature distribution. 

Taguchi's tactics for robust design are to specify a criterion 
for assessing the effect of the noise factors and to estimate it 
by experimentation. While noise factors are difficult or imprac- 
tical to control in production or in use, for purposes of the ex- 
periment (to learn about the effect of the noise factors), the noise 
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factors need to be controlled during the experiment. The s- 
expected loss L{ a criterion for assessing the effect of the 
noise factors at a particular combination of control factor levels 
■^control* can be defined for a general loss function /(-) [29]: 

^ (^control) = i ^(^(•r C ontiol^noi 8 e))*/(X 

noise) ^*noiie* (2) 

Notation 

/( •) loss function 

^control. *noi«) response at combination of control factor 

levels (x MIIIroI> x^J 
£(*comroi) ^-expected loss at 
/( -) joint pdf of noise factors. 

The objective of robust design is to find a product/process 
design, x^*^, with minimum ^-expected loss. An appropriate 
loss function for reliability is discussed in section 3.3. 

3.1 Product Array Designs 

Taguchi [26] proposed estimating the loss (2) via ex- 
perimentation according to specialized experimental plans (pro- 
duct arrays). A product array consists of 2 plans or arrays, one 
for the control factors (control array) and the other for the noise 
factors (noise array). The product-array design is so named 
because all the noise-factor combinations specified by the noise 
array are run with every combination of the control factors 
specified by the control array. 

As an example, consider an experiment to improve the 
reliability of a night-vision goggle tube sub-assembly [19]. The 
tube sub-assembly is insulated by a combination of protective 
coatings which degrades over time and over exposure to humidi- 
ty & temperature. One goal of the experiment was to make the 
tube insulation reliability robust to the handling of tube when 
the coatings are applied. The noise factor, A/, had 2 levels (the 
tube is handled or not). Ten control factors (A - J) were chosen, 
from many, across a 17-step manufacturing process and include 
factors related to the tube packaging such as type of coating, 
primer, and electrical-connection configuration. Each control 
factor was also studied at 2 levels according to a 12-run Plackett- 
Burman design, eg, 2 types of primer coating (factor B) and 
2 types of lead coating material (factor H) were used. Therefore, 
the product array consisted of a 12-run Plackett-Burman design 
for the control factor array and a simple noise factor array (1 
factor at 2 levels) as displayed in table 6. At each of the 24 
control/noise factor combinations, one tube was manufactured 
and then life tested under a high-temperature cycling and humidi- 
ty regimen. The tubes were inspected for failure every 2 days 
for 20 days. The lifetime data in table 6 make some assump- 
tions because it is unclear from [19] whether 16 of the tubes 
failed between days 18-20 or whether they were still function- 
ing at 20 days. For illustration, they are treated here as still 
functioning, resulting in right-censored data; 6 of the tubes also 
failed before the inspection #1 at day 2, yielding left-censored 
data. The product-array design used in this experiment is a 
fractional-factorial design which would otherwise have required 
2 = 2048 control and noise factor settings. 



TABLE 6 

Product-Array Design & Lifetime Data 
[Goggle Experiment] 



Control Array Noise Array, N 



A 


B 


C 


D 


£ 


F G 


H 


/ 


/ 


1 


2 


1 


1 


1 




1 


1 1 


1 


1 


1 


(0,2) 


(20,oo) 


1 


1 


1 


1 


2 


2 2 


2 


2 


2 


(20,oo) 


(20,oo) 


1 


1 


2 


2 


1 


I 1 


2 


2 


2 


(0,2) 


(0.2) 


1 


2 


1 


2 


I 


2 2 


1 


1 


2 


(7,9) 


(20, oo) 


1 


2 


2 


1 


2 


1 2 


1 


2 


1 


(20,«) 


(20,oo) 


I 


2 


2 


2 


2 


2 1 


2 




1 


(20,oo) 


(7.9) 


2 


1 


2 


2 


1 


2 2 


I 


2 


1 


(20,oo) 


(20,oo) 


2 


1 


2 


1 


2 


2 1 


1 


1 


2 


(20,00) 


(0.2) 


2 


1 


1 


2 


2 


1 2 


2 


1 


1 


(0.2) 


(0,2) 


2 


2 


2 


1 


1 


1 2 


2 


1 


2 


(20,oo) 


(20,oo) 


2 


2 


1 


2 


2 


1 1 


1 


2 


2 


(20,oo) 


(20,oo) 


2 


2 




1 


1 


2 1 


2 


2 




(20,oo) 


(20,00) 



3.2 Analysis of Product Array Design Data 

For analyzing the product-array data, Taguchi [26] pro- 
posed estimating the j-expected loss ^(x^,) (2) for each 

^control 

specified by the control array from the data obtained by 
varying the noise factors according to the noise array and then 
modeling the estimated losses in terms of the control factors. 
That is, he proposed constructing signal-to-noise ratios from 
the noise-array data and analyzing them by common methods 
for designed experiments such as analysis of variance. Ref [21] 
referred to Taguchi's approach as the loss-model approach. 
Alternatively, [29] proposed modeling the response Kdirecdy 
in terms of both the control & noise factors and then evaluating 
the loss using the estimated response model. Their rationale for 
the latter (response-model [21]) approach was that it would be 
more likely to find a simple model for the response than one 
for the much more complicated estimated s-expected loss. Ex- 
amples in [21, 29] gave evidence for preferring the response- 
model approach because it provides additional insight. 

For reliability applications, the response-model approach 
is reasonable because the same parametric regression models 
in section 2 can be used. An even more compelling reason 
against using the loss-model approach is that the signal-to-noise 
ratios cannot handle censored data. The product-array data allow 
model (1) to be fit consisting of all Control main effects (with 
possibly some Control x Control interactions), all 
Control x Noise interactions and all Noise main effects (with 
possibly some Noise x Noise interactions), where Control & 
Noise denote specific control & noise factors, respectively. The 
covariate values have the form: 

X ~~ ^comtoIi *controlX noise* X noisc ), 

for the main effect values, 

(•^control » noise) 
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and corresponding Control x Noise interaction values, 

'coniroix noise* 




C1 
C2 



N 

Figure 1. Example Response-Functions Showing Opportuni- 
ty for Robustness 

The Control x Noise interactions are important in achiev- 
ing robust reliability because the fact that the ^-expected loss 
(2) changes for various control-factor combinations means that 
these interactions must exist. Figure 1 displays a simplified rela- 
tionship between a response Kand 1 control-factor C (at 2 levels) 
and 1 noise-factor N (varied over an interval). The vertical lines 
in figure i show that the effect of the noise factor on the response 
Yis substantially smaller at control-factor level 1 (CI) as the 
noise factor varies and hence is more robust at that level. Con- 
sequently, robust design exploits the existence of such interac- 
tions between control- and noise- factors. 

3.3 Choosing Parameter Values (Factor Levels) 

Once estimates for the response-model effects have been 
obtained, the important control-factor levels need to be recom- 
mended. For a simple model with few noise factors, recom- 
mended levels might be apparent from inspecting the model 
direcdy; ie> by observing what the important effects. Shoemaker 
ex al [21] gave an example, but for complicated models, their 
approach might be tedious. 

An alternative is to specify some meaningful criterion or 
s-expected loss (2) and use the identified response model (1) 
to evaluate them, assuming some distribution for the noise fac- 
tors. In practice, where it is difficult to specify a distribution 
for the noise factors, the criterion can be evaluated over a sam- 
ple of noise combinations as specified by a noise array. The 
same noise array from the experiment need not be used, 
however. For example, instead of using a fractional-factorial 
design, the s-expected loss could be evaluated over a full- 
factorial design. Similarly, the s-expected loss could be 
evaluated for all possible settings of the control factors, not just 
those specified by the control array. 



For achieving robust reliability, therefore, reliability should 
depend as little as possible on the noise factors. Also high 
average reliability is required. Ref [6] considered criteria bas- 
ed on the linear part (x T -fi) of model (1), viz, the mean 
log(lifetime) for the lognormal regression model. This paper 
assesses reliability in terms of the probability of exceeding a 
certain time 7*, such as a warranty period. Using model (1), 
this survival probability is; 

Sf{(tog(D ~ x T >0)/o}. (3) 
Notation 

Sf{} an appropriate survivor function 

(li 'control » •'control X noise* ■'noise) • 

Thus, using (3) as the loss function in (2) means that the s- 
expected loss here is the ^-expected reliability. For a given 
control-array combination, the j-expected reliability can be 
estimated by evaluating the reliabilities (3) over all the noise- 
array combinations and then taking their mean. These evaluated 
reliabilities (over the noise array) represent a sample of 
reliabilities which can be summarized in other meaningful ways; 
eg, the standard deviation which measures the variability of the 
reliabilities as the noise factors vary. By taking a worst-case 
approach, the minimum of the sample can also be compared. 
Based on these criteria, control-array combinations with large 
mean, large minimum, and small standard deviation are 
desirable. Analysis of the goggle experiment in section 7 il- 
lustrates the use of these criteria. 

3.4 Usual Taguchi Robust Reliability-Design 

Most of the robust design applications involve a continuous 
quality characteristic Y with ideal or target value Y 0 and use 
the squared error loss function (Y- Y 0 ) 2 in (2). Provided fac- 
tors can be found which can adjust the mean quality 
characteristic to target, then robust design amounts to minimiz- 
ing variability of the quality characteristic caused by the varia- 
tion of the noise factors. Situations in which reliability is defined 
in terms of the quality characteristic exceeding a threshold fall 
into this category. For example, Taguchi often gives the exam- 
ple of an electronic system in which the parameters of the com- 
ponents drift with age. In this example, the control factors levels 
are nominal values for the component parameters and the noise 
factors levels are deviations in the component parameters from 
the nominal values that reflect the drifting. No lifetesting is re- 
quired here, which is a distinct advantage of analyzing the quali- 
ty characteristic direcdy. This paper focuses on situations which 
do not have such continuous quality characteristics, however, 
so that lifetesting must be done. 

4. ANALYSIS METHODS FOR CENSORED DATA 

For model (1), when all the lifetimes are observed (com- 
plete sample), the analysis is straightforward using ML estima- 
tion [12]. The MLE for (/?, o) are found by maximizing the 
following likelihoods: 
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For the Weibull regression model, 

£03. o) = Yl OM-exp[[(>v - xf.0]/a] 

i 

- expKa - x?.0)/o]l 
For the lognormal regression model, 

i 

The estimates are compared with their standard errors [12] to 
identify the important effects. The problems arising from cen- 
sored data are discussed next. 

One method which unwisely continues to be used treats 
the right-censoring times as actual failure times and then 
analyzes them by common methods for complete data. (For 
interval-censored data, an interval endpoint or midpoint might 
be used.) Although simple, ignoring the censoring can lead to 
wrong decisions because the unobserved failure times and right- 
censoring times can differ greatly, depending on the particular 
factor-level combination. A simulation [8] showed that this 
method can perform quite poorly by missing some important 
effects and mis-identifying spurious effects. 

ML estimation easily handles both failure and censored 
(right or left) data. The MLE for (jS, a) are found by maximiz- 
ing the following likelihoods: 
For the Weibull regression model, 

£(0, o) = Yl (l/*)-«p{[Ov - xT.0]/a\ 

i € FAIL 

-exp[(y ( -jr/.0)/a]} ■ exp{exp[-( Vj - x[.&)td\}. 

f€CENS 

(4) 

For the lognormal regression model, 

£<0, o) = Yl (l/a)-gaud[(y f - xJ-$)/o\ 

/€ FAIL 

• Yl g aufc [U ~ *r (5) 
ieCENS 

Standard errors for the MLE can be obtained [12]. Com- 
mercially available software perform these computations, such 
as SURVIVAL, the SYSTAT survival analysis module [24], 
or the LIFEREG procedure in SAS [20]. For a censored datum, 
its contribution to £ is simply the probability of being censored. 
Similarly for an interval-censored datum (a,b); its contribution 
to £ is the probability of failing between times a & b. 

One problem with ML estimation for censored data is that 
the MLE might not exist, ie, at least one parameter estimate 
is infinite, so that testing for important effects cannot be done 



by comparing the MLE with their standard errors. Ref [22] gave 
necessary & sufficient conditions for the existence of MLE for 
model (1). In the reliability context, [7] concluded that 
estimability problems tend to occur for the designs in sections 
2 & 3 where the fitted model has nearly the same number of 
parameters as number of observations. 

The estimability problem of the ML approach motivated 
the Bayes approach [10]. The Bayes approach is reasonable 
because important factor effects might be anticipated to be large 
but not infinite. By using appropriate prior distributions, well- 
behaved posterior distributions can be calculated whose finite 
modes or medians can be used as estimates of the effects. Also, 
posterior distributions allow the importance of factorial effects 
to be assessed without using the asymptotic approximations 
needed by the ML method; a simple way to assess whether a 
factor effect is important is to see whether its marginal posterior 
does not contain zero. Ref [10] considered the lognormal regres- 
sion model (1) and used the reasonable conjugate prior [18]; 
the prior pdf is: 

p(A a) = a-*.exp(-(j3 - 0o) r '4><0 - j3 0 )/(2a 2 )] 

.a" (,,0+I) .exp[-^ 0 .^/(2a 2 )]. (6) 

The posterior is proportional to the product of the likelihood 
(5) and the prior (6) and is relatively simple to obtain numerical- 
ly using recent advances in Bayes computing. An appropriate 
choice of the prior parameters (v 0 , si, (3$, Aq) gives a very dif- 
fuse or essentially non-informative prior which means that the 
results reflect the information provided by the data. The sen- 
sitivity of the results based on the current data can be assessed 
by trying more informative priors [10]. The Bayes approach 
is illustrated in section 7 in the analysis of the goggle experiment. 



5. ANALYSIS OF THE FLUORESCENT-LAMP 
EXPERIMENT 

Consider the fluorescent-lamp experiment in section 2. The 
experiment studied 5 factors (A - E) using the design in table 
2. Besides the 5 main effects (A- E) t the experimenters thought 
that the A xB interaction might be important. Using ML estima- 
tion, a lognormal regression model was fit using the lifetime 
data in table 2. Table 7 gives the MLE and ^-significance levels 
(p values) for the 5 main effects ( A - E) and the A xB interac- 
tion (the intercept is denoted by INT). Based on these results, 
the main effects D, B, E, A are important (small p values) in 
the order given. Therefore, only 4 of the 5 factors are impor- 
tant, with A being only marginally important. Since the AxB 
interaction is not important, the signs of the important main ef- 
fects can be used to recommend factor levels. If the effect is 
negative, then level 1 (whose covariate value is -1) improves 
reliability, viz* the mean log(lifetime). Similarly, for a positive 
effect, level 2 (whose covariate value is + 1) improves reliabili- 
ty. Thus the results from table 7 suggest that reliability can be 
improved at factor levels A^BiD^E^ where the subscript in- 
dicates the recommended level. 
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Table 7 

MLE and P Values of lognormal Regression Model 
[Fluorescent Lamp Experiment] 



Effect 


MLE 


P Value 


INT 


2.939 


0.000 


A 


-0.117 


0.059 


B 


0.201 


0.001 


AB 


-0.049 


0.430 


C 


0.051 


0.408 


D 


-0.273 


0.000 


E 


0.153 


0.015 


a 


1.590 


0.000 



6. ANALYSIS OF THE THERMOSTAT EXPERIMENT 

Consider the thermostat experiment in section 2 which 
studied 1 1 factors {A - K) . Taking the ML approach, a lognor- 
mal regression model with 1 1 factor main effects {A - K) was 
fit using the lifetime data (table 3); the results are in table 8 
(model 1). (At most 12 effects can be fit simultaneously because 
of the design run size of 12, so that besides an intercept only 
the main effects could be considered in an initial analysis.) Based 
on table 8, 9 of the factors appear to be important. One reason 
for so many apparently important factors was pointed out by 
Bullington et al [3]: Each group of 10 units was produced at 
the same time so that the variability among the 10 units tends 
to be smaller than if they had been produced at different times. 
This reduced variability, which is used in the statistical testing, 
is one possible explanation for many ^-significant effects. Some 
additional analysis based on [8, 9] which account for the pro- 
perties of the 12-run Plackett-Burman design used in this ex- 
periment, suggests the presence of an ExH interaction. Fur- 
ther evidence of an interaction is seen in table 8 (Model 1) by 
noting that the MLE for all factors except E & H have nearly 
the same magnitude. Consequently, model 2 was fit in which 
the factor B main effect (the least important from model 1) was 
dropped and replaced by the ExH interaction. The results in 
table 8 (model 2) indicate that only the E, H, ExH effects are 
important. The original analysis [3] found only E & H impor- 
tant and recommended E X H\. The same recommendation is ob- 
tained using the table 8 results. While the original analysis did 
not account for the possibility of an interaction, the same recom- 
mendations result because the two factors have a synergistic 
relationship, Wz, both increase reliability as they are changed 
from their level 2 to their level 1 . In summary, the analysis sug- 
gests that reliability gains can be made by using a grain size 
(factor E) of 0.008 inches with a heat treatment (factor H) of 
45 minutes. 



7. ANALYSIS OF THE GOGGLE EXPERIMENT 

In the robust reliability experiment to improve night-vision 
goggles in section 3, there were 10 control-factors (A - J) and 
1 noise-factor (N). Taking the response-model approach, a 



lognormal regression model (1) can be fit which consists of an 
intercept (INT), 10 Control main effects (A - K), 1 Noise main 
effect {N), and 10 Control x Noise interactions (AxN - KxN), 
where Control & Noise denote control & noise factors, respec- 
tively. Using the product-array data in table 6, the MLE for 
this model do not exist. Consequently, the Bayes approach [101 
was taken. Using a relatively diffuse prior, table 9 gives the 
central 0.95 & 0.99 intervals of the marginal posteriors for each 
effect. The important effects appearing in bold face are those 
whose central 0.95 intervals do not contain zero; in fact, the 
0.99 intervals for all of these except the IxN interaction do not 
contain zero. Based on these results, main effects, B, D - H y 
and interactions, AxN t CxN, ExN, IxN, JxN, are important. 



TABLE 8 

MLE and P Values of Lognormal Regression Models 
[Thermostat Experiment] 





Model I 






Model 2 




Effect 


MLE 


P Value 


Effect 


MLE 


P Value 


A 


-0.312 


.0001 


A 


-0.091 


.3890 


B 


0.221 


.0024 


EH 


0.663 


.0024 


C 


-0.319 


.0001 


C 


-0.098 


.3474 


D 


0.285 


.0001 


D 


0.064 


.5174 


E 


-1.023 


.0001 


E 


-1.023 


.0001 


F 


0.231 


.0016 


F 


0.010 


.9219 


G 


-0.390 


.0001 


G 


-0.169 


.1075 


H 


-0.557 


.0001 


H 


-0.557 


.0001 


I 


-0.332 


.0001 


I 


-0.112 


.2872 


J 


-0.277 


.0001 


J 


-0.056 


.5958 


K 


-0.352 


.0001 


K 


-0.131 


.2149 



For this experiment, the relationship between the response 
and the control & noise factors is too complicated to make con- 
trol factor-level recommendations simply by inspecting the 
model. Consequently, the criteria in section 3 based on the 
reliability (3) distribution (mean, standard deviation, and 
minimum) can be evaluated over the 2 levels of the 1 noise- 
factor for each of the possible combinations of control factors 
(1024 = 2 10 ) and then ranked appropriately (out of 1024, with 
I being the best). At each control factor setting, the reliability 
(3) was evaluated using a time T = 100 with (/3, a) being sampl- 
ed 5000 times from the posterior distribution at each noise 
factor-level which yielded a total of 10 4 reliabilities. Thus, the 
uncertainty in (0, a) as described by the posterior distribution 
is easily incorporated by the Bayes approach, which is an at- 
tractive advantage. Table 10 presents the 15 best control-factor 
combinations according to the mean criterion, and the standard- 
deviation & minimum criteria. Based on these results, a good 
choice of factor levels is: 

Afi 2 C x D x E 2 E 1 G x H 2 I?J2 (row #\ of table 10). 

Its minimum criterion is the best and its standard-deviation 
criterion is also small which indicate that these levels are robust 
to the noise factor. Contrast this with the settings in row #1 
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of table 10 which has a relatively high mean but is clearly not 
robust to the noise factor. 



TABLE 9 

Posterior Quantiles Using Diffuse Prior 
[Goggle Experiment] 



Quantile 



Effect 


.005 


.025 


.975 


.995 


INT 


2.60 


2.64 


2.99 


3.07 




-0. 15 


-0.08 


0. 32 


0.40 


B 


0.74 


0.80 


1.15 


1.24 


C 


-0.23 


4). 18 


0.23 


0.28 


D 


-0.79 


-0.71 


-0.32 


-0.24 


E 


0.31 


0.37 


0.77 


0.84 


F 


0.04 


0.10 


0.52 


0.60 


G 


-0.60 


-0.52 


-0.12 


-0.07 


H 


0.27 


0.32 


0.70 


0.82 


I 


-0.35 


-0.27 


0.12 


0.20 


J 


-025 


-0.18 


0.23 


0.30 


N 


-0.11 


-0.05 


0.29 


0.35 


AN 


-0.68 


-0.62 


-0.21 


-0.16 


BN 


-0.12 


-0.05 


029 


0.35 


CN 


-0 ,82 


-0.72 


-0.34 


-0.27 


DN 


-0.22 


-0.17 


0.24 


0.31 


EN 


-0.58 


-0.50 


-0.07 


-0.01 


FN 


-0.09 


-0.03 


0.39 


0.45 


GN 


-0.46 


-0.37 


0,03 


0.09 


HN 


-0.35 


-0.24 


0.16 


0.23 


IN 


-0.53 


-0.44 


-0.04 


0.03 


JN 


-0.79 


-0.71 


-0.33 


-0.27 


a 


0.02 


0.02 


0.07 


0.11 
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2 
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11 
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2 
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2 


1 


1 2 


2 
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11 




2 
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2 
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2 2 
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2 
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15 
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DOE/Opt: A System for Design of Experiments, 
Response Surface Modeling, and Optimization 
Using Process and Device Simulation 

Duane S. Boning, Member, IEEE, and P. K. Mozumder, Member, IEEE 



Abstract — Rapid modeling and optimization of manufacturing 
processes, devices, and circuits are required to support modern 
integrated circuit technology development and yield improve- 
ment. We have prototyped and applied an integrated system, 
called DOE/Opt, for performing Design of Experiments (DOE), 
Response Surface Modeling (RSM), and Optimization (Opt). The 
system to be modeled and optimized can be either physical 
or simulation based. Within the DOE/Opt system, coupling to 
external simulation or experimental tools is achieved via an 
embedded extension language based on Tel. The external problem 
then appears to DOE/Opt as a model with user defined inputs and 
outputs. DOE/Opt is used to generate splits for experiments, to 
dynamically build and evaluate regression models from experi- 
mental runs, and to perform nonlinear constrained optimizations 
using either regression models or embedded executions. The 
intermediate regression modeling can appreciably accelerate the 
optimization task when simulation or physical experiments are 
expensive. The primary application of DOE/Opt has been to 
process optimization using coupled process and device simulation. 
DOE/Opt has also been applied to process and device simulator 
tuning, and to aid in device characterization. Such a DOE/Opt 
system is expected to augment the use of TCAD tools and to 
utilize data collected by CIM systems in support of process 
synthesis. We have demonstrated the application of the system 
to process parameter determination, simulator tuning, process 
control modeling, and statistical process optimization. We are 
extending the system to more fkilly support emerging device 
design and process synthesis methodologies. 



I. Introduction 

SEMICONDUCTOR process and device design can sub- 
stantially benefit from the use of simulation and modeling 
to reduce the cost and time required to develop new or extend 
existing technology [11. Technology development, however, 
requires substantially more than a fundamental simulation 
capability: tools and methods to assist in exploration of design 
trade-offs and to optimize a design are becoming increasingly 
important. 

Work on frameworks for the integration of technology CAD 
(TCAD) tools has focused primarily on underlying representa- 
tions for the wafer and process, and on the interfaces between 

Manuscript received August 1, 1993; revised December 7, 1993. This work 
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Microelectronics Technology Office under Contract F33615-88-C-5448. 
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P. K. Mozumber is with the Semiconductor Process and Design Center, 
Texas Instruments. Inc., P.O. Box 655012, Dallas, TX 75235 USA. 
IEEE Log Number 9216485. 



simulation tools and these representations l2]-f4]. Systems 
which achieve integration and help automate the execution 
of tools have also been demonstrated [5]-[7], Less common 
are tools or systems which directly address how one uses such 
integrated simulation capability to solve design problems [8]. 
Alvarez explored statistical modeling methods for the design 
and optimization of processes and devices [9], In that work, 
formal (Box— Behnken) design of experiments were used to 
construct response surface models. A grid search method was 
used to explore the design space and to identify feasible 
regions given goals and constraints on multiple performances 
and performance sensitivities. Other design exploration and 
optimization techniques have been investigated elsewhere, 
particularly for circuit yield optimization [10] — [12]. 

This paper makes contributions in two areas. First, we 
describe approaches for the use of experimental design, re- 
gression or direct simulation modeling, and optimization that 
we have found to be important in process/device design, 
simulator tuning, coupling to process control, and design for 
manufacturability. Second, we describe a solution to issues 
in the implementation of a general purpose optimization tool 
suitable for end TCAD users including graphical interface, user 
programmability via an extension language, and interfaces to 
existing simulation tools. DOE/Opt is a "task level" tool which 
complements existing TCAD Framework research on "data 
level" representations, and builds on "tool level" simulation 
capability [13]. Two audiences should thus benefit from this 
paper: 1 ) designers who will increasingly demand the types of 
simulation, modeling, and optimization capability represented 
by DOE/Opt, and 2) implementors who will find it necessary 
to implement similar high level utility programs in the future. 

In Section II, the overall architecture of the prototype 
DOE/Opt system is introduced. Each of the key conceptual 
components of the system is discussed in succeeding sections. 
The encapsulation and use of simulation or other models, and 
the construction of analytic response surfaces is discussed in 
Section III. The integration of formal experimental design 
techniques with 1 ) simulation or experimental execution, and 
2) design exploration and model construction are described in 
Section IV. The integration and use of nonlinear optimization 
is discussed in Section V. Four examples are presented in 
Section VI to illustrate the use of the DOE/Opt system. These 
include process/device performance optimization, simulator 
model tuning, process control recipe generation, and manufac- 
turability optimization. Conclusions are drawn in Section VII. 
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Fig. 1. The architecture of the DOE/ Opt system. 

IL System Architecture 

The DOE/Opt system architecture is summarized in Fig. 1 . 
There are three key computational components to the system: 
I) design of experiments (DOE); 2) response surface model 
generation (RSM); and 3) nonlinear constrained optimization 
(Opt). The DOE/Opt executable uses the Tel [14] extension 
language for application and end-user programmability. On 
top of this sits a Motif compatible graphical user inter- 
face (GUI) implemented in Tk [15] to manipulate models. 
Through the GUI, the end-user directs the execution of model 
blocks, constructs response surface models, and defines and 
directs optimizations. Optimizations are performed by a sec- 
ond executable program which implements multiple objective 
nonlinear constrained optimization on top of the NPSOL [16] 
package. DOE/Opt consists of approximately 10000 lines of 
Tcl/Tk code, and about 1 000 lines of C code (to interface to 
NPSOL). 

A typical DOE/Opt work flow used to explore or optimize 
design trade-offs is shown in Fig. 2. The first task is to 
automate, via a block script, the simulation and data extrac- 
tion flow. Second, one chooses an appropriate experimental 
design to systematically vary block inputs, and performs the 
corresponding simulations. Third, one constructs and evaluates 
response surface models for the outputs. Additional design 
points may need to be simulated and different modeling 
strategies employed to achieve an adequate model fit. Once 
achieved, target values and constraints on the outputs are 
defined, and optimization using the response surface models 
(or using the body directly) are performed. One or more 
candidate optima are verified using the full simulation, which 
may indicate the need for additional simulation and model 
building in the region of interest. 

In the following sections, the key components in the 
DOE/Opt architecture will be described in more detail. In 
Section VI, variations on the usage flow of Fig. 2 will be 
discussed for a number of optimization scenarios. 

III. Model Representation 

Experimental design and optimization fundamentally re- 
quire a way to express the response of the system under 



Automate simulation flow 

- write block body script 

- define block inputs and outputs 

\ 

Explore Design Space 

- choose experimental design 

- perform simulations 

1 Augment design 

Build RSMs 
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for better model 
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Perform Optimization ! 

- define ta^sar^r^ints SSSSwmwte 

- perform multiple optimizations ^ , 



Verify Results 
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of interest 



Results Adequate? 



Fig. 2. Row diagram of typical DOE/Opt system use. 



mos.block 




Fig. 3. The conceptual organization of a DOE/Opt block. 

study. Conceptually, a generic representation of a model of that 
system is needed, independent of whether analytic, numeric, 
or experimental execution of the model produces the response. 
In this way, many of the same experimental design, execution, 
and optimization mechanisms can be applied to real systems, 
to extensive numerical simulations, or to fast analytic models. 

In DOE/Opt, a block is the embodiment of some model, 
and maps inputs to outputs. A block has several tabular and 
script components as pictured in Fig. 3. The key computational 
component is the block body. The block body is often a script 
that chains together process simulation, device simulation, 
and output data extraction. In other cases, a response surface 
model (or RSM block) contains a block body which is an 
analytic polynomial function to calculate response and gradient 
information. Similar to the block body is the block context, 
which is another script that is run only once when the block 
is first loaded. The context is useful for defining utility 
procedures that will be called repetitively within the block 
body. 
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Fig. 4. The DOE/Opi user interface. Pictured is a block lor the calculation and optimization of key 256 Mb DRAM pass transistor parameters, including 
display of inputs, outputs (targets and constraints), a portion of the simulation script in Ihe block body, and optimization results. 



The block input table defines the possible inputs to the 
block. For each input, the name, default value, units, minimum, 
and maximum values may be specified. A toggle for each 
input indicates whether or not that input will vary in an 
experimental design or optimization. The block coefficient 
table is very much like the input table, except that it defines 
the values of any internal model coefficients inside the block. 
For example, automatically generated response surface models 
have coefficient values determined by the regression fit to run 
data. The block output table defines the possible outputs that 
the block is able to compute. Each output has a name, and may 
have an associated "RSM" giving the filename of an external 
response surface model that may be used in place of the block 
body to calculate the output. Desired transformations of output 
values may also be specified. For optimization, additional 
information about each output may be specified in the output 
table: each output may have a target value, and/or may have 
upper and lower constraint values. Finally, each row in the 



run and optimization tables presents specific input values and 
resulting outputs (or the initial starting point and resulting 
optima in the case of the optimization table). 

A. Concrete Representation 

We represent each block as a separate UNIX file. The block 
file contains information that may be loaded into the DOE/Opt 
user interface where it can be viewed, edited, executed, and 
saved back to a UNIX file. The DOE/Opt graphical user 
interface for a typical block is shown in Fig. 4. This block file 
is also directly executable and manipulatable from the UNIX 
command line; this is crucial for the use and manipulation 
of blocks in other scripts, programs or systems external to 
DOE/Opt. For example, a block which chains process and 
device simulation to calculate threshold voltage might be 
executed from the UNIX command line, where input values 
are specified as command line arguments and the result is 
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returned by the block: 

unix> vt .block -ox-time 100 

-channel-length 0.5 

0.675 

Alternatively (and more commonly), the same execution might 
be performed directly from the user interface. In this case, one 
would define an experiment via the "Run/Design" pulldown 
menu which would generate new rows in the run table. 
The user might optionally override the inputs, coefficient, or 
run table entries for the ox_time and channel .length 
variables, and then execute the run table entries via the 
"Execute" pulldown menu. In this case, the block body would 
be run and the results displayed in the output columns of the 
run table. 

The block may be manipulated in other ways, and an 
entire set of runs and optimizations (as defined per some 
experimental design) can be invoked, either from the graphical 
user interface or from the UNIX command line. For example, 
to change the target value for the threshold voltage output, and 
then to perform the optimization specified inside the block, one 
could use the unix command line commands: 

unix> vt. block set output vt -target 0.7 
unix> vt .block optimize.rsm. 

In the graphical user interface, one would simply point to the 
entry widget in the target column for the vt output row 
to change the value, and then initiate the optimization via the 
"Optimize using RSMV option in the "Optimize" pulldown 
menu. This performs the optimization using response surface 
models for the desired outputs (vt in this case), and fills in 
the optimization table with the optimal operating points. In the 
UNIX invocation, the resulting optimization table will reside 
within the block file at the completion of the optimization. 
To support long simulation or optimization problems, it is 
critical that the graphical user interface not be required; a 
UNIX command line approach, for example, enables overnight 
runs without a graphical display. 

B, Encapsulation of Simulation Flows 

We have made two specific uses of the generic block inter- 
face. One is to represent automatically constructed response 
surface models as described later in Section IV-C. The other 
is to automate simulation and data extraction [17]. To support 
the user definition of simulation flows, the following four 
issues are important. First, a high-level interpreted scripting 
language is essential for dynamic generation and use of 
scripts. We use the Tel (Tool Command Language [14]) 
extension language for all block body scripts, thus ensuring 
that full language mechanisms (e.g., procedure definitions, 
variable manipulations, mathematical expression evaluation, 
and control Mow) are available as needed. Second, the script 
must be well connected to the surrounding block, as that 
is where inputs, coefficients, and outputs are specified. We 
accomplish this by making all input and coefficient values 



specified via the GUI (as shown in Fig. 4) available as 
Tel variables in the block body script, and by returning 
requested outputs back from the script to the GUI. Third, 
"template" versions of simulation input files are important. The 
parameterization of simulation decks, along with parameter 
substitution to produce actual simulation input files, is common 
practice [9]. We extend this so that full expression evaluation, 
and indeed Tel procedure invocation, is also possible within 
the template. For example, a SUPREM3 [18] deck might be 
parameterized as 

Initialize silicon (100) boron=lel5 
Thickness=2 

Diffusion t ime=#[$ox_time*60 . 0 ] # 
temp=1000 dry02 

Implant arsenic dose=lel5 energy=$= 

{energy} 

so as to convert time in hours to time in minutes for use 
in the simulation. An interactive utility program to convert 
conventional input files to DOE/Opt templates has also been 
implemented, where selected ranges of text are turned into 
input variables and the appropriate input item is dynamically 
added to the block definition in the DOE/Opt GUI. Finally, 
we have found that interactive control over block execution 
is important. A "monitor" pop-up can be made to appear for 
selected simulations within a block body or for the block body 
execution as a whole. The monitor shows the current execu- 
tion status and permits termination of problem executions or 
optimizations. 

IV. Design of Experiments 

The value of statistically based experimental designs (the 
matrix of runs generated by specific combinations of input 
values) has been well established [9], [19]. The DOE/Opt 
system strongly encourages the use of designed simulation ex- 
periments where the input parameters are varied in a carefully 
structured pattern so as to maximize the information that can 
be extracted from the resulting simulations. In our work, we 
have achieved tight integration between specification of mod- 
els (via blocks), generation of experimental designs, automatic 
execution of those designs, and generation of response surface 
models from execution results. 

A. Experimental Design Suite 

A variety of experimental designs are incorporated into 
DOE/Opt [19]. A nominal design generates a single run with 
nominal (default) values for all selected block inputs. A 
center point design is similar, but uses the midpoint value 
(between the maximum and minimum limits) for all selected 
inputs. Note that the range for input values may be specified 
directly, as deltas from the default value (e.g., "+5"), or as 
percentages from the default value (e.g., "-10%"). The axial 
design includes center points, and points at the minimum and 
maximum of each input while holding other inputs to their 
nominal values. The Box-Wilson design includes the center 
point, and axial points and corner points distributed around an 
n -dimensional sphere circumscribing the cuboid defined by the 
input minima and maxima [20). We have found two variants 
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of the conventional Box-Wilson design to be useful. The Box- 
Wilson (inscribed) modifies the design by distributing axial 
and corner points on an n-dimensional sphere inscribed within 
the cuboid [21]. This has the benefit that no experimental 
points will be outside the input ranges, but results in models 
that are not well defined in the corners of cuboid. The Box- 
Wilson on a cube stretches the comer points to lie on the 
cuboid comers to achieve more complete coverage of the input 
parameter space. The full-factorial design produces a uniform 
grid with user-specified density covering the input parameter 
space. Other conventional designs, including Box-Behnken 
and half fractional (as well as the half fraction complement 
for augmenting a screening design), are provided [19]. Finally, 
Latin hypercube sampling (LHS) [22J provides an orthogonal 
array that randomly samples the entire design space broken 
down into r n equal-probability regions (where r is the number 
of runs, and n is the number of input variables). LHS can be 
looked upon as a stratified Monte Carlo sampling where the 
pairwise correlations can be minimized to a small value (which 
is essential for uncorrelated parameter estimates) or else set to 
a desired value [23]. LHS is especially useful in exploring the 
interior of the parameter space, and for limiting the experiment 
to a fixed (user specified) number of runs. All designs are 
generated algorithmically, with the exception of Box-Behnken 
designs which use table lookups. Latin hypercube samples 
are generated via an interface to LHS software from Sandia 
National Laboratories [24]. 

B. Integration with Block Execution 

When an experimental design is selected, the DOE/Opt sys- 
tem examines the block to determine which inputs have been 
selected to vary, and to determine the input defaults and limits. 
Based on the design a run table is generated (or added to) with 
empty slots for each of the outputs to be calculated. The run 
table can then be executed in total or incrementally to fill in 
missing elements. The run table may be edited after its creation 
and before executing the block, so that the user has control 
over the simulation run value, and can delete or add new run 
table rows to eliminate or perform other simulations. After 
execution, the run table can be examined using interactive 
data visualization methods. In the current implementation we 
interface to the xgobi package, an interactive dynamic graphics 
program which is effective in the exploration and visualization 
of multivariate data [25]. Finally, response surface models can 
be automatically generated from the run table. 



C. Response Surface Models 

An important part of the DOE/Opt system is the generation 
and use of response surface models to inexpensively "mimic" 
the more complex workings of a simulation or experimental 
block. The following four elements of response surface model 
construction have been found to be important. First, a variety 
of polynomial base models (linear, quadratic with and without 
cross terms, cubic, etc.) must be provided. We have found it 
useful to link the regression model type with the experimental 
design used (e.g., one builds a "Box-Wilson" response surface 
model from run data generated via the Box- Wilson design). In 



addition, scaling of the inputs to the range [-1, 1] is used to 
increase model construction robustness. The covariance of the 
estimates, a metric of model stability, is dependent on the input 
design matrix and the lack of model fit. Choice of the input 
design matrix is critical to determining the model coefficients, 
and minimizing the covariance between the estimates of the 
model coefficients. Scaling the inputs minimizes the corre- 
lation between the estimates of the coefficients of the model 
[26], [27]. Additionally, transformation of outputs may be used 
(e.g., log, exp y square, square root, inverse) to aid model fitting 
and later optimization. Second, response surface models must 
be capable of generating both response and response gradient 
information in order to reduce evaluations required during 
optimization. The automatically generated response surface 
models have as their output the basic response, followed by the 
Jacobian with respect to the input parameters, followed by the 
Jacobian with respect to the model coefficients (to aid in model 
tuning if desired). Code to calculate all of these analytically is 
part of the generated response surface model. Third, weighted 
regression is important: a common sequence is to perform a 
broad experimental design, build a model, optimize to some 
region within that model, then refine the experimental design 
near that region. In such cases, it is useful to weight the 
second set of simulations more than the first to increase model 
accuracy in the region where the model will be most used. 
The response surface models are generated using weighted 
least squares regression, currently via an interface to the Splus 
statistical package [28], [29]. Finally, tight integration between 
the generated model and DOE/Opt is necessary. This includes 
the ability to perform executions and optimizations using either 
a full simulation or using generated response surface model 
blocks. It is also important that response surface models be 
generated on the fly: a compile, link, and re-execution of the 
application is unacceptable. For this reason, we again elected 
to generate interpretive, Tcl-based regression models. 

V, Optimization 
A key role of the DOE/Opt system is to provide an easy to 
use optimization capability. Real design and manufacturability 
optimization problems often have the following features. First, 
the optimization problems are often nonlinear, and multiple 
targets or objectives must be reconciled simultaneously. Sec- 
ond, numerous constraints apply in the optimization problem; 
these constraints may also be nonlinear in nature. Many design 
and optimization problems focus on choices of continuous 
parameter values. While discrete value or integer optimization 
is also often desirable (e.g., to make a choice in whether 
to include a design feature or not, or to decide between 
two alternative treatments), we have focused entirely on the 
continuous parameter optimization problem. 

We have constructed a layered optimization capability as 
shown in Fig. 5. The fundamental optimization capability is 
provided by the NPSOL package [16]. On top of this a set 
of Tel bindings provide interpretive access to the FORTRAN 
NPSOL routines. In addition to enabling one to dynamically 
"call into" the optimizer, the optimizer also "calls out" to 
the Tel language layer for the evaluation of objective and 
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Fig. 5. Layering of the Tct/NPSOL interface. The higher components make 
use of the lower level blocks. The shaded components are described in this 
document. 
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Fig. 6. Detail of the calling sequence in the Tel to NPSOL interface. 



constraint functions, as illustrated in Fig. 6. In this way it 
becomes possible to define an objective procedure and call 
the optimizer (which in turn evaluates the objective function 
several times), all dynamically and at run time. Layered on 
top of the basic access to NPSOL via Tel, a generic weighted 
multiple target, multiple constraint optimization capability 
has been implemented. This layer is responsible for folding 
together and making use of all existing model and model 
gradient information (such as that provided by the DOE/Opt 
generated response surface models) to reduce optimization 
solution times. Through this interface, DOE/Opt is able to 
generate models and optimization problems, and solve them 
via NPSOL. 

Several key issues have been identified and addressed in 
order to make the tool more appropriate for use in semiconduc- 
tor process and device design. First, in practice optimization 
problems can become highly nonlinear and complex. No 
guarantees on global optimality are made. Instead, an approach 
to easily support the optimization from numerous starting 
points is used. In this case, we use the same design of 
experiments methods outlined in Section IV to produce well 
distributed coverage of the parameter space in searching for 
potential optima. The Latin hypercube sample has proven 
particularly useful in generating optimization starting points. 

Second, we have found that a selection of overall objective 
functions needs to be provided. The current system includes 
the following choices; weighted sum of squares, weighted 



sum of normalized squares, weighted sum of absolute er- 
rors, weighted cumulative errors, and weighted normalized 
cumulative error. While some of these are not necessarily 
well-behaved (e.g., sum of absolute errors), they often map 
well onto the conceptual optimization problem at hand and 
work satisfactorily. 

Third, we have found that it usually makes sense for 
the optimization to use transformed response surface model 
information rather than the "pure" values returned by models 
constructed under transformations. For example, a log trans- 
form selected when building a response surface model might 
fit 10 14 , 10 15 , and 10 16 as a straight line internally. However, 
the model will still return exponentially large values when 
evaluated. For optimization, on the other hand, the linear 
underlying model is more effectively used. The ability to 
perform model regressions and optimizations in transformed 
space is important. 



VI. Examples 

Four examples demonstrate the application of the DOE/Opt 
system. These include process/device design, simulator tuning, 
process control recipe generation, and statistical process/device 
design, 

A. Process/ Device Optimization 

A very common problem in coupled process and device 
design is to choose process parameters so as to achieve desired 
device performance. DOE/Opt has been applied to study the 
selection of channel doping in the design of a 256Mb DRAM 
pass transistor. The choice of energies and doses for threshold 
adjust and punchthrough implants (vtudose, vtjenergy,ptjdose t 
ptjenergy) are the input parameters under study. The per- 
formance parameters are leakage current (Ueakage), peak 
channel doping (Njnax) y and charge voltage (V .charge, or V D 
at which I D > I .critical, with V G = 3.5 V, V s = 2.5 V, and 
Vb = -1-5 V), which is a measure of the rate at which the 
storage capacitor can be charged. 

A simulation script executes the SUPREM3 process sim- 
ulation and MEDICI [30] device simulation sequence given 
a set of process parameter values; part of this script is 
shown in Fig. 4. A Box-Wilson on a cube experimental 
design was first run over a relatively broad input parameter 
range, and quadratic models constructed. The V. charge and 
Njnax models fit the data well. The good N jnax fit is 
expected, since the maximum doping concentration in a single 
implant is directly proportional to implant dose and inversely 
proportional to implant straggle (which is only a weak function 
of implant energy). The V .charge fit is also not surprising, 
as there is an inverse dependence of V .charge on V G - V T 
for small Vds, and V T is roughly linearly dependent on 
implant dose; the resulting inverse dependence is weak enough 
for our conditions that it can be well approximated by a 
quadratic response surface model. However, the Ueakage 
model required the use of logarithmic transforms, and still 
achieved a relatively poor fit due to highly nonlinear behavior 
when punchthrough occurs. By weighting the simulation runs 
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prior to regression, an adequate fit near the I .critical parameter 
was achieved. 

Based on these models, an optimization was performed 
with a high Vxharge output selected as the target (where we 
minimized the squared error from the target), with constraints 
on Ueakage and Njnax as indicated in Fig. 4). The resulting 
trade-off between Ueakage and Vxharge is clearly seen in 
Fig. 7, which also shows the original design points and 
candidate optima (possible optimal points suggested by the 
nonlinear optimization using the response surface models) 
for the threshold adjust parameters. Additional iterations of 
simulation, model building, and optimization as in Fig. 2 could 
be used to obtain and verify optima to a desired level of 
accuracy. An interesting conclusion based on study of the two 
implants is that both implants share functionality in setting 
the threshold and protecting against punchthrough. Future 
work using DOE/Opt will study the necessity of two channel 
implants as compared to a single channel adjust implant. We 
also note that the feasible region appears to be narrow; the 
need for manufacturability analysis as described in Section 
VI-D is clearly indicated. 

5. Simulator Tuning 

DOE/Opt has been applied to simulator model tuning. In 
this case, the goal was to determine ionization coefficients 
in an internally enhanced version of PISCES (31] to match 
measurements of substrate current. The model is of the form 
a = a°° exp [-(E CTlt /E n ) n ]. PISCES simulation decks were 
developed which calculate substrate current (7 5 ) at drain volt- 
ages V D = 5, 6, and 7 volts, and gate voltage V G = 1, 2, 3, and 
4 volts. The input parameters to the script are the ionization 
model coefficients a°°, E C[it , and exponent n. A Box-Wilson 
experimental design varying these parameters was performed, 
and regression models of Is(V G ,V D ) as a function of the 
ionization model coefficients were constructed. Optimization 
was performed using the response surface models to find 
ionization model coefficients that achieved the best fit to 
the experimental data (minimum sum of squared error). The 
result was approximately 15% average error across several 
V D values. A comparison between the substrate currents from 
measurements, simulation before calibration, and simulation 
after calibration is made in Fig. 8. 

Two observations on the use of DOE/Opt in this problem 
bear special mention. First, we found that DOE/Opt is well 
suited to optimization and fitting to a small and mixed number 
of targets and constraints. In its current form, it is not ideal 
for "parameter extraction" where one desires to fit to large sets 
of data. Work to improve DOE/Opfs capability by interfacing 
to the MACH system is underway [32]. Second, we find that 
the most difficult task facing a user of DOE/Opt is writing the 
data extraction routines used within a script body. We found 
it necessary to develop and make available to users utility 
routines that extract data (e.g., current- voltage curves) from 
process and device simulation output files. To support higher 
level tools such as DOE/Opt, it is important that the textual and 
binary output formats of simulation tools be well documented 
and available. 
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Fig. 7. Exploration of 256 Mb pass transistor channel implant design. 
The trade-off between leakage and charging voltage is shown in (a). The 
corresponding threshold adjust design points are shown in (b). The original 
Box-Wilson experimental design points are shown as circles (feasible points 
satisfying the design constraints are shown as open circles), and the candidate 
points from optimization are shown as pluses. 

C. Process Control Recipe Generation 

DOE/Opt has been used to generate optimal recipes for 
a plasma enhanced chemical vapor deposition of silicon ni- 
tride (PECVD nitride) process run on the Applied Materials 
Precision 5000 reactor (AMT5000). As a part of the Micro- 
electronics Manufacturing Science and Technology (MMST) 
program our aim was to model (physically and empirically), 
optimize, and adaptively control the process to target [33], 
[34]. We have used DOE/Opt to generate the initial recipe* 
for the process. Using the models of the PECVD nitride 
process, the optimal equipment settings for the AMT5000 were 
generated to obtain the desired outputs from the process. 
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Fig. 8. Comparison of experiments substrate current, modeled substrate 
current using calibrated model coefficients, and modeled current using default 
(pre-calibrated) model coefficients. 

The DOE/Opt setup for recipe generation is shown in Fig. 
9. The equipment controllables are specified via the input 
table. The primary equipment controls on the AMT 5000 for 
the PECVD nitride by the nitrogen-silane-ammonia process 
are the N 2 flow (7V 2 ). SiH 4 flow (SiH 4 ), NH 3 flow (AW 3 ), 
pressure (/V), RF power (F) t and electrode gap (Gap). The 
ranges are chosen by the process engineer based on his 
understanding of the limitations of the equipment and validity 
of the process models. The default values are the centers of 
the hyperbox defined by the ranges. The quality characteristics 
of interest include the film deposition rate {rate), index of 
refraction (ri), stress {stress), and thickness nonuniformity 
(nu). These responses, the corresponding desired values, and 
the specification limits are indicated in the output table. The 
block body contains the Tel script used to encapsulate the pre- 
existing PECVD nitride model program (invoked using the C 
executable pecvd). The body script handles the conversion 
to and from that expected by the pecvd program, 

A Box-Wilson on a cube experimental design was used to 
create full quadratic response surface models. The run table 
comprised 77 rows, each of which corresponded to a point 
in the input design space. For each of the rows the model 
was evaluated and all four output values were obtained by 
executing the body. Once the run table is filled, regression 
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Fig. 9. DOE/Opt setup for recipe generation of the PECVD nitride process 
on the AMT5000. 

is used to generate full quadratic response surface models 
(consisting of 28 coefficients) for each of the outputs in terms 
of all six inputs. The generated RSM block files are linked to 
the optimization problem via the "RSM" column in the output 
table. 

A weighted sum of squares optimization is performed to 
determine the optimal recipe for the process. The weights 
for the outputs are chosen as the inverse of the standard 
errors of the models (see [34] for rationale). For the PECVD 
nitride process rate, stress, and ri are targets, whereas nu acts 
as a constraint (the specific target and constraint values are 
shown in Fig. 9). The objective function is specified to be a 
weighted sum of squares of the difference between the model 
prediction and the target. In one instance of recipe generation, 
all of the inputs were allowed to vary for the optimization, 
optimal points were found and verified, and the controller was 
initiated with the generated initial recipe; the controller was 
subsequently able to adjust the recipe to response to equipment 
shifts and drifts [34], A second instance of recipe generation 
is described here, in which the process engineer wanted to 
generate an optimal recipe where the Gap was not changed 
from its value of 0.42 (either for the initial recipe or during 
control). Minimal changes had to be made to the DOE/Opt 
block to generate the initial recipe under the new constraint. 
The "Vary" button corresponding to the Gap in the input table 
was deselected (to toggle the "vary" value for the gap in Fig. 
9). A new set of starting points, without the variable Gap, was 
then generated for the optimizer. To minimize the probability 
of getting a solution at a local minimum, multiple starting 
points were used for the optimization: the nominal point and a 
10 point Latin hypercube design was used. Corresponding to 
each of the starting points an optimal point was generated. Two 
distinct optima were found where the objective function values 
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Fig. 10. Verification of optima] recipe for the PECVD nitride process on the 
AMT5000 with Gap held constant. 
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Fig. 11. Schematic of manufacturability simulation. 



are relatively small: for N 2 = 1800. SiH 4 = 160, NH 3 = 49 
(with error 0.025) and N 2 = 2200 (with error 0.028). Both of 
the optima are verified using the PECVD nitride model, and 
the results are shown in Fig. 10. Run #2 was chosen as the 
best recipe since it had a closer to target value of ri by trading 
off a small amount on stress. The controller was activated 
using the initial recipe generated by DOE/Opt. The controller 
is presently in routine use [34]. 

D. Design for Manufacturability 

DOE/Opt was used with a statistical process simulator 
to study design for manufacturability. Several criteria for 
optimizing process recipes under random process variations 
have been described in the literature [35]-[39J. The concepts 
of parametric yield, standard deviation and "signal to noise 
ratio" have been implemented using DOE/Opt and used for 
process optimization. We present an example of gate oxi- 
dation process optimization under simulated manufacturing 
variations. Using this example, we compare three different 
design for manufacturability strategies. 

The semi-analytic statistical process and device simulator 
FABRICS [40], [41] is used to emulate the effects of process 
and material parameter disturbances on the resulting device 
structure. Given the nominal settings for a process step and 
the standard deviations for the process disturbances, we seek 
to determine the mean shift in the disturbances (or equivalently 
new process parameter nominal values) that will result in a 
design that is optimal under some manufacturability criterion. 

We encapsulated FABRICS into DOE/Opt via a body Tel 
script which generates FABRICS input files, runs the simula- 
tor, and parses the resulting output file. For each DOE/Opt run, 
FABRICS is executed in the statistical mode to generate 100 
samples, from which the statistical response parameters are 
calculated. As shown in Fig. 11, our approach is to construct 
response surface models for each of several responses as 
a function of the inputs and disturbances. In the oxidation 
example below, a Box-Wilson on a cube design was executed 
and full quadratic models for each output were generated using 
least squares regression. We found that all of the models had 
excellent goodness of fits. 

Gate Oxidation Problem Formulation: Fig. 12 shows the 
DOE/Opt problem formulation for optimizing a gate oxidation 
process. The set up is analogous to that described in the 



previous section. The input table consists of six variables, 
There are two process parameters: time of gate oxidation 
(time^ox), and temperature of gate oxidation (temp.gox). 
Associated with each of the process parameters are the mean 
and the standard deviation for the disturbances, denoted by the 
prefixes rnd. and sd.. The means and the standard deviations 
for the disturbances, which are assumed to be independent and 
normally distributed, can be determined by tuning FABRICS 
to the fabrication line for which the process is to be optimized 
|42]. In using DOE/Opt, we have explicitly specified the means 
and standard deviations of the disturbances corresponding to 
time^gox and temp^ox, and do not vary the standard deviations 
(i.e., they are fixed at their tuned values). The output table 
consists of five parameters. These correspond to the mean of 
the gate oxide thickness {mean Jhsiol gate), standard deviation 
of the gate oxide thickness (sigjhsio2gate\ signal-to-noise 
ratio of the gate oxide thickness (s2nJhsio2gate) defined as 
the ratio of the mean to the standard deviation, mean of 
Monte Carlo parametric yield for the gate oxide thickness 
(ymJhsio2gate\ and the standard deviation of Monte Carlo 
parametric yield for the gate oxide thickness iysdJhsio2gate). 
The standard deviation is used to determine the confidence 
interval on the yield based on the Monte Carlo sampling. Yield 
calculation is based on whether the resulting responses fall 
within a region of acceptability; the limits for the acceptability 
of the gate oxide thickness are specified using the DOE/Opt 
coefficients table. 

Gate Oxidation Optimization Results; Three different ob- 
jective functions were used as manufacturability criteria for 
the gate oxidation optimization. The aim was to observe if 
one particular metric was better than the others for solving 
the problem of statistical recipe generation. Following are the 
objective functions used as the test cases: 
Case 1: Min sigjhsio2gate (std. dev.) 
Case 2: Max s2nJhsio2gate (signal to noise) 
Case 3: Max ymjhsio2gate (Monte Carlo yield) 

s.t. ysdjhsio2gate < 5% 
All cases: s.t. 4. 1 < mean Jhsio2 gate < 4.3 
The minimization and maximization were carried out using 
a small and a large target value, respectively. For each case, 
optimizations (smallest sum of square errors from defined 
targets) were carried out from multiple starting points. The 
verification runs of the best optima for each case are shown 
in Fig. 12. 

It was found that for each of the three cases the optima 
for all the starting points resulted in very nearly the same 
values for the outputs. Moreover, the final values of the output 
parameters are almost identical irrespective of the objective 
function. On closer observation it is also found that the 
optimal process inputs, mdJime^ox and mdJemp.gox, are 
identical for the first two cases, i.e., where sigJhsio2gate and 
$2nJhsio2gate are used as objective functions. The optimal 
process inputs when yield was used as an objective function 
are slightly different. 

These results are shown pictorially in Fig. 13. The nearly 
identical answers in all three cases was due to the stringent 
constraint on mean Jhsio2 gate. All three metrics improve 
with smaller values for the standard deviation of the gate 
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Inputs 



vary 


name 


default 


unit* 


min 


max 


0 




1.600 


see*le3 


1.000 


2.200 


0 




1.250 


K*le3 


1.160 


1.360 


1 




0 


sec+le3 


-0.1 


0.1 


0 


sd-time^ox 


0 


sec*le2 


0 


0.1 


1 


md_t*mp_gojt 


0 


K*le3 


-0.1 


0.1 


0 


sd_temp_gox 


0 


K*le2 


0 


0.1 



Coefficients 



tune 


name 


value 


min 


max 


1 


Lim_thiio2gate 


4.2 


3.7 


4.7 



Outputs 



■elect 


name 


target 


units 


min 


max 


weight 


rim 


Constraint 


mean_thsio2gate 


4.2 


lc-8m 


4.1 


4.3 




mean thsio2gate.rsm 


TWget 


sig_thsio2gatc 


le-2 


lo-8m 


1.0e-2 


1 




sigL.thsio2gate.rsm 




■2n_thsio2g*tc 


100 




1 


100 


1.0 


s2iL.thsio2gate.rBm 




Ym.thsio2gate 


100 


% 


50 


100 


1.0 


vm-thsio2gate.rsm 


Constraint 


ysd_thsu>2gate 


0 


% 


0 


6 




ysd-thsio2gate.rsm 



Opt Method Comparisons 



case 


rnd, 
timc-gox 


md. 
tcmp_gox 


mean_ 
thsio2gate 


th*io2gate 


»2n_ 
thsio2gate 


ym_ 

thsio2gate 


y«d- 

thsio2gate 


1 
2 
3 


-0.1 
-0.1 
-0.0165 


0.O1O67 
0.010960 
-0.0041596 


4.300021075 
4.300007041 
4.182424063 


0.3063008676 
0.3063571703 
0.3191762864 


14.14464137 
14.13267310 
13.25157406 


84.34461003 
84.38151207 
87.1727678 


3.187603261 
3.187613642 
3.477796522 



Fig. 12. DOE/Opt setup for statistical process optimization of gate oxide thickness using FABRICS, and verified simulation results for three opti- 
mization strategies. 

oxide (sigjhsio2gate). Yield is maximized if in addition to 
a decreasing variability the design is "centered" — the mean 
value is chosen so that most of the distribution is within 
the constraints of acceptability. Similarly the signal-to-noise 
is maximized as the mean value increases without adversely 
affecting the standard deviation, or the standard deviation 
decreases without a significant drop in the mean value. For 
the case where the mean value is tightly constrained (towards 
the center of the window in this case) all three metrics 
are optimized when sigJhsio2gate is minimized, subject to 
the constraints on mean Jhsiol gate. The constraint on the 
mean ensured that both s2njhsio2gate and ymjhsio2gate 
are strongly dominated by sigJhsio2gate. This is the case 
even when the sigJhsio2gate and mean Jhsio2 gate are not 
independent. The difference in the results between the yield 
maximization and the remaining metrics can be attributed to 
the fact that the distribution of the gate oxide thickness results 
in larger yield loss in the tail regions when signal to noise or 
standard deviation alone are optimized. 

The dependence of s2nJhsio2gate and ymJhsio2gate on 
meanjhsio2gate was explored by repeating the same optimiza- 
tions without the constraints on the mean Jhsio2 gate. Results 
(shown schematically as dashed curves in Fig. 13) indicated 
that both s2nJhsio2gate and ymJhsio2gate were strong func- 
tions of both sigjfisio2gate and mean Jhsio2 gate. Moreover, 
the sensitivities of the yield and signal to noise ratios to the 
mean and the standard deviation were significantly different. 
Maximizing ymjhsio2gate produced the same results with 
and without the constraint on mean Jhsio2 gate. On the other 




Fig. 13. Manufactur ability optimization for 1) minimum standard deviation 
(constrained mean in solid lines, unconstrained in dashed line), 2) maximum 
signal to noise, and 3) maximum yield. 

hand, the sigjhsio2gate nunirnization and s2nJhsio2gate max- 
imizations produced identical results which were significantly 
different from the corresponding results with the constraints 
on mean Jhsio2 gate. The optimizations were dominated by 
decreasing the variability of the gate oxide thickness and the 
resulting process conditions produced extremely low yields, 
i.e., the optimizer was able to push the process conditions 
near the edge of the acceptability region. 
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This example demonstrates that DOE/Opt can be used to 
assist in statistical process design, an integral part of design for 
manufacturability. Several different objective functions can be 
defined and the results of using different statistical criteria for 
optimization can be explored. We have been able to derive and 
calculate multi-dimensional yield by encapsulating a statistical 
process simulator, FABRICS. Finally, we have been able to 
determine the process conditions that maximize the value of 
the yield using the optimizer NPSOL in DOE/Opt. 

VII. Conclusion 

A system has been implemented to provide design of exper- 
iments, regression modeling, and optimization capability for 
use in Technology CAD. The system emphasizes the uniform 
treatment of model evaluation across numerical simulation 
and response surface modeling. In this fashion, effective 
optimization layered on top of existing simulation capability is 
possible. We have demonstrated the application of the system 
to process parameter determination, simulator tuning, process 
control modeling, and statistical process optimization. Future 
work is needed to extend the system to more fully support 
emerging device design and process synthesis methodologies 
[43]. 
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